{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trip Purpose Estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:516: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint8 = np.dtype([(\"qint8\", np.int8, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:517: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_quint8 = np.dtype([(\"quint8\", np.uint8, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:518: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint16 = np.dtype([(\"qint16\", np.int16, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:519: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_quint16 = np.dtype([(\"quint16\", np.uint16, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:520: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint32 = np.dtype([(\"qint32\", np.int32, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorflow/python/framework/dtypes.py:525: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  np_resource = np.dtype([(\"resource\", np.ubyte, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:541: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint8 = np.dtype([(\"qint8\", np.int8, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:542: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_quint8 = np.dtype([(\"quint8\", np.uint8, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:543: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint16 = np.dtype([(\"qint16\", np.int16, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:544: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_quint16 = np.dtype([(\"quint16\", np.uint16, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:545: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint32 = np.dtype([(\"qint32\", np.int32, 1)])\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:550: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  np_resource = np.dtype([(\"resource\", np.ubyte, 1)])\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import shapely\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import collections\n",
    "import sklearn\n",
    "import sklearn.preprocessing\n",
    "import sklearn.model_selection\n",
    "import sklearn.ensemble\n",
    "import rfpimp\n",
    "from sklearn.feature_selection import SelectFromModel\n",
    "import random\n",
    "import datetime\n",
    "\n",
    "# oversampling\n",
    "from imblearn.over_sampling import SMOTENC\n",
    "from imblearn.under_sampling import RandomUnderSampler\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data should contain the nearby POI information from the Google Place API. \n",
    "\n",
    "`Preprocessing is needed`\n",
    "\n",
    "Consider using the POI types from the paper by Emagun et al."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_16 = pd.read_csv('../../Data/model_inputs/gdf_2016_X.csv')\n",
    "y_16 = pd.read_csv('../../Data/model_inputs/gdf_2016_y.csv')\n",
    "X_17 = pd.read_csv('../../Data/model_inputs/gdf_2017_X.csv')\n",
    "y_17 = pd.read_csv('../../Data/model_inputs/gdf_2017_y.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_all = X_17#pd.concat([X_16, X_17],axis=0)\n",
    "y_all = y_17#pd.concat([y_16, y_17],axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The encoding of **purpose** and **mode** is as follows. Note that the trips of home and work have been removed.\n",
    "\n",
    "- purpose: {'leisure': 0, 'food_drink': 1, 'shops': 2, 'pick_up_drop_off': 3, 'education': 4, 'health': 5}\n",
    "- mode: {'walking': 0, 'public_transport': 1, 'car': 2, 'cycling': 3, 'public_transport, car': 4, 'other': 5}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2    9906\n",
       "0    5913\n",
       "3    3514\n",
       "1    3380\n",
       "4     175\n",
       "5      60\n",
       "Name: mode, dtype: int64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_all['mode'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['id_trip', 'mode', 'duration', 'distance_m', 'weekday', 'precip',\n",
       "       'temp', 'morning', 'midday', 'afternoon', 'evening', 'midnight',\n",
       "       'startx', 'starty', 'endx', 'endy'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_all.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The meaning of X variables is shown below.\n",
    "\n",
    "| Variable | Meaning | Unit |\n",
    "|:-------- | :-------- | :---- |\n",
    "| id_trip | Unique ID | NA |\n",
    "| mode | Trip mode | {'walking': 0, 'public_transport': 1, 'car': 2, 'cycling': 3, 'public_transport, car': 4, 'other': 5} |\n",
    "| duration | Trip duration | Second |\n",
    "| distance_m | Euclidean trip distance calculated from nodes along a Shapely LineString object | Metre |\n",
    "| weekday | Whether this trip occurs in a weekday (Monday to Friday) | {1:True, 0:False} |\n",
    "| precip | Total precipitation within a given hour on a given date | mm |\n",
    "| temp | Average Temperature within a given hour on a given date | Celcius |\n",
    "| morning | Whether this trip occurs in the morning time 0600-1059 | {1:True, 0:False} |\n",
    "| afternoon | Whether this trip occurs in the afternoon time 1100-1359 | {1:True, 0:False} |\n",
    "| evening | Whether this trip occurs in the evening time 1700-2159 | {1:True, 0:False} |\n",
    "| midnight | Whether this trip occurs in the midnight time 2200-0559 | {1:True, 0:False} |\n",
    "| startx | The x coordinate of the starting point in the TBC reference system | Metre |\n",
    "| starty | The y coordinate of the starting point in the TBC reference system | Metre |\n",
    "| endx | The x coordinate of the ending point in the TBC reference system | Metre |\n",
    "| endy | The y coordinate of the ending point in the TBC reference system | Metre |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace False and True with 0 and 1, respectively\n",
    "X_all = X_all.replace({True: 1, False: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id_trip</th>\n",
       "      <th>mode</th>\n",
       "      <th>duration</th>\n",
       "      <th>distance_m</th>\n",
       "      <th>weekday</th>\n",
       "      <th>precip</th>\n",
       "      <th>temp</th>\n",
       "      <th>morning</th>\n",
       "      <th>midday</th>\n",
       "      <th>afternoon</th>\n",
       "      <th>evening</th>\n",
       "      <th>midnight</th>\n",
       "      <th>startx</th>\n",
       "      <th>starty</th>\n",
       "      <th>endx</th>\n",
       "      <th>endy</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>150744</td>\n",
       "      <td>2</td>\n",
       "      <td>862</td>\n",
       "      <td>9935.922336</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>16.910884</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7.624322e+06</td>\n",
       "      <td>1.247673e+06</td>\n",
       "      <td>7.631864e+06</td>\n",
       "      <td>1.250415e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>35763</td>\n",
       "      <td>2</td>\n",
       "      <td>1873</td>\n",
       "      <td>6832.113937</td>\n",
       "      <td>1</td>\n",
       "      <td>9.972328e-07</td>\n",
       "      <td>18.007062</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7.641919e+06</td>\n",
       "      <td>1.236661e+06</td>\n",
       "      <td>7.641941e+06</td>\n",
       "      <td>1.236732e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>323826</td>\n",
       "      <td>2</td>\n",
       "      <td>1165</td>\n",
       "      <td>12233.968564</td>\n",
       "      <td>1</td>\n",
       "      <td>9.972328e-07</td>\n",
       "      <td>18.007062</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7.625699e+06</td>\n",
       "      <td>1.248427e+06</td>\n",
       "      <td>7.616671e+06</td>\n",
       "      <td>1.252917e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>375668</td>\n",
       "      <td>2</td>\n",
       "      <td>820</td>\n",
       "      <td>3612.405991</td>\n",
       "      <td>1</td>\n",
       "      <td>1.495849e-06</td>\n",
       "      <td>20.940647</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7.630948e+06</td>\n",
       "      <td>1.253394e+06</td>\n",
       "      <td>7.630941e+06</td>\n",
       "      <td>1.253445e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>212877</td>\n",
       "      <td>2</td>\n",
       "      <td>1316</td>\n",
       "      <td>26428.385191</td>\n",
       "      <td>1</td>\n",
       "      <td>1.495849e-06</td>\n",
       "      <td>23.092603</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7.631318e+06</td>\n",
       "      <td>1.252962e+06</td>\n",
       "      <td>7.646416e+06</td>\n",
       "      <td>1.264378e+06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   id_trip  mode  duration    distance_m  weekday        precip       temp  \\\n",
       "0   150744     2       862   9935.922336        1  0.000000e+00  16.910884   \n",
       "1    35763     2      1873   6832.113937        1  9.972328e-07  18.007062   \n",
       "2   323826     2      1165  12233.968564        1  9.972328e-07  18.007062   \n",
       "3   375668     2       820   3612.405991        1  1.495849e-06  20.940647   \n",
       "4   212877     2      1316  26428.385191        1  1.495849e-06  23.092603   \n",
       "\n",
       "   morning  midday  afternoon  evening  midnight        startx        starty  \\\n",
       "0        1       0          0        0         0  7.624322e+06  1.247673e+06   \n",
       "1        0       1          0        0         0  7.641919e+06  1.236661e+06   \n",
       "2        0       1          0        0         0  7.625699e+06  1.248427e+06   \n",
       "3        0       1          0        0         0  7.630948e+06  1.253394e+06   \n",
       "4        0       0          1        0         0  7.631318e+06  1.252962e+06   \n",
       "\n",
       "           endx          endy  \n",
       "0  7.631864e+06  1.250415e+06  \n",
       "1  7.641941e+06  1.236732e+06  \n",
       "2  7.616671e+06  1.252917e+06  \n",
       "3  7.630941e+06  1.253445e+06  \n",
       "4  7.646416e+06  1.264378e+06  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_all.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['purpose'], dtype='object')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_all.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The size of each trip purpose:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0    8049\n",
       "2    7723\n",
       "1    2714\n",
       "4    2175\n",
       "3    1379\n",
       "5     908\n",
       "Name: purpose, dtype: int64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('The size of each trip purpose:')\n",
    "y_all['purpose'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variable settings and algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " There are four settings of variables\n",
    " - Basic (without nearby trip purposes or Google nearby Places)\n",
    " - Basic and nearby trips\n",
    " - Basic and Google POIs\n",
    " - Baisc and nearby trips and Google POIs \n",
    " \n",
    " \n",
    " There are two algorithms:\n",
    " - Random Forest\n",
    " - Nested Logit\n",
    " \n",
    " \n",
    " Based on exhaustive combinations of settings and algorithms, eight models will be trained.\n",
    " \n",
    " In this notebook, we will compare four random forest models. Each model will be tuned using GridSearchCV.\n",
    " \n",
    " The basic variables include the following:\n",
    " - 'mode'\n",
    " - 'duration'\n",
    " - 'distance_m'\n",
    " - 'weekday'\n",
    " - 'precip'\n",
    " - 'temp'\n",
    " - 'morning' \n",
    " - 'midday' \n",
    " - 'afternoon'\n",
    " - 'evening'\n",
    " - 'midnight'\n",
    " \n",
    " \n",
    " There is no need to normalise continuous variables before using RF."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A function that runs a standard random forest model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_rf(X_train, X_test, y_train, y_test, n_estimators=100, max_depth=None):\n",
    "    \"\"\"\n",
    "        Run a random forest classification model\n",
    "    \"\"\"\n",
    "    clf = sklearn.ensemble.RandomForestClassifier(n_estimators=n_estimators, max_depth=max_depth)\n",
    "    clf.fit(X_train, y_train)\n",
    "    score = clf.score(X_test,y_test)\n",
    "    preds = clf.predict(X_test)\n",
    "    #print(\"Random Forest Classifcation accuracy:\", score)\n",
    "    return clf, score, preds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A function to compute the proportions of nearby trip purposes in the training data (TK)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_nearby_purposes(y_train, X_train, X_test, buffer_size_m=200):\n",
    "    \"\"\"\n",
    "        Augment the training and testing data with the proportions of nearby trip purposes in the training data.\n",
    "        :param y_train: (dataframe) the trip purpose of the training data, which is used to compute proportions of nearby trip purposes. Should contain a column of 'purpose'\n",
    "        :param X_train: (dataframe) the predictor variables of the training data. Should contain three columns of 'id', endx' and 'endy'.\n",
    "        :param X_test: (dataframe) the predictor variables of the testing data. Should contain two columns of 'id', 'endx' and 'endy'.\n",
    "        :param buffer_size_m: (int) the size of the buffer in metres around each trip end.\n",
    "        :return: (list) [X_train_aug, X_test_aug]. X_train_aug and X_test_aug are the X_train and X_test augmented with extra columns, respectively. The new columns are ['p_nearby_0', 'p_nearby_1', 'p_nearby_2', 'p_nearby_3', 'p_nearby_4', 'p_nearby_5']\n",
    "        Note: when selecting the trips close to a trip in the training data, remember to exclude this trip itself using the 'id' column.\n",
    "        Note: EPSG:3347 == Canada Lambert Projection\n",
    "    \"\"\"\n",
    "    ## merge train and test X back together\n",
    "    all_X = pd.concat([X_train, X_test],axis=0)\n",
    "    \n",
    "    ## initialise a geo-dataframe of all the X values and create a buffer of 200 m around each trip (for each row)\n",
    "    all_X['geometry'] = all_X.apply(lambda row: shapely.geometry.Point(row['endx'],row['endy']),axis=1)\n",
    "    geo_X = gpd.GeoDataFrame(all_X,crs=\"EPSG:3347\")\n",
    "    \n",
    "    # calculate 200 m buffers around each trip end point\n",
    "    geo_X['buffers'] = geo_X['geometry'].apply(lambda row: row.buffer(buffer_size_m))\n",
    "    \n",
    "    # get a geo-dataframe with only the Trip ID,trip end buffer and trip purpose\n",
    "    only_buffers = geo_X[['id_trip','buffers']]\n",
    "    only_buffers = gpd.GeoDataFrame(only_buffers.rename(columns={'buffers':'geometry'}), crs=\"EPSG:3347\")\n",
    "    \n",
    "    \n",
    "    ## intialise a geo-dataframe of only the training data\n",
    "    # Note: only the purposes of trips of the training dataset are known so this is why we need to create a seperate geoDataFrame of only the training data\n",
    "    all_train = pd.concat([X_train,y_train],axis=1)\n",
    "    all_train.reset_index(drop=True,inplace=True)\n",
    "    all_train['geometry'] = all_train.apply(lambda row: shapely.geometry.Point(row['endx'],row['endy']),axis=1)\n",
    "    geo_train = gpd.GeoDataFrame(all_train,crs=\"EPSG:3347\")\n",
    "\n",
    "    ## perform the spatial join between a buffer of each trip and trip end points from exclusively the training data.\n",
    "    # Note: here we are extracting the purpose of trips where the trip end points from the training geo-dataframe (type==Point) intersect with a buffer of all the trip end points (type==Polygon) \n",
    "    joined_data = gpd.sjoin(only_buffers, geo_train, op='intersects', how='left')\n",
    "\n",
    "    ## drop duplicates (as each trip will fall within a buffer of itself)\n",
    "    to_drop = joined_data[['id_trip_left','id_trip_right']].apply(lambda row: True \\\n",
    "                                                                if row['id_trip_left'] == row['id_trip_right']\\\n",
    "                                                                else False, axis=1)\n",
    "    joined_data = joined_data[~to_drop]\n",
    "    \n",
    "    ## Compute the proportion of nearby purposes types for each trip/row\n",
    "    grouped_data = joined_data.groupby(['id_trip_left', 'purpose']).agg({'mode':'count'})\n",
    "    grouped_data = grouped_data.unstack().apply(lambda row: row/row.sum(),axis=1)['mode'].reset_index()\n",
    "    grouped_data = grouped_data.fillna(0.0)\n",
    "    \n",
    "    ## rename columns\n",
    "    new_columns = ['id_trip','p_nearby_0','p_nearby_1','p_nearby_2','p_nearby_3','p_nearby_4',\\\n",
    "                       'p_nearby_5']\n",
    "    grouped_data.columns = new_columns\n",
    "    \n",
    "    ## Merge p_nearby columns back to X_train and X_test\n",
    "    all_X = all_X.merge(grouped_data, on='id_trip', how='left')\n",
    "    all_X = all_X.fillna(0.0)\n",
    "    X_train_aug = X_train.merge(all_X[new_columns], on='id_trip')\n",
    "    X_test_aug = X_test.merge(all_X[new_columns], on='id_trip')\n",
    "    return X_train_aug, X_test_aug"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic RF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "basic_variables = ['mode', 'duration', 'distance_m', 'weekday', 'precip',\n",
    "       'temp', 'morning', 'midday', 'afternoon', 'evening', 'midnight']\n",
    "\n",
    "# needed for nearby purpose computation \n",
    "temporary_variables = ['id_trip','endx','endy']\n",
    "\n",
    "# extract all variables\n",
    "X_all_basic_variables = X_all[basic_variables + temporary_variables]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(22948, 14)\n",
      "(22948, 1)\n",
      "(22948,)\n"
     ]
    }
   ],
   "source": [
    "print(X_all_basic_variables.shape)\n",
    "print(y_all.shape)\n",
    "print(np.ravel(y_all).shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# split training testing data\n",
    "rd_state = 42\n",
    "\n",
    "X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(\n",
    "    X_all_basic_variables, y_all, stratify=y_all, random_state=rd_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# compute nearby purposes\n",
    "X_train, X_test = compute_nearby_purposes(y_train, X_train, X_test, 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preserve Trip IDs of training and testing data\n",
    "X_train_IDs = X_train['id_trip']\n",
    "X_test_IDs = X_test['id_trip']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# remove columns which are not needed\n",
    "X_train = X_train.drop(columns=temporary_variables, axis=1)\n",
    "X_test = X_test.drop(columns=temporary_variables, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# flatten y train data\n",
    "y_train = np.ravel(y_train)\n",
    "y_test = np.ravel(y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the random forest algorithm on this dataset\n",
    "rf_clf, rf_accuracy_rest, rf_pred_test = run_rf(X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy on training data:  0.9999418976236129\n",
      "Accuracy on testing data:  0.5438382429841381\n"
     ]
    }
   ],
   "source": [
    "# print results\n",
    "print(\"Accuracy on training data: \", rf_clf.score(X_train, y_train))\n",
    "print(\"Accuracy on testing data: \", rf_accuracy_rest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prediciton accuracy on the training and testing data is 1.0 and 0.511, respectively, indicating that the model is overfitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.54      0.67      0.60      2012\n",
      "           1       0.35      0.11      0.17       678\n",
      "           2       0.54      0.63      0.58      1931\n",
      "           3       0.45      0.15      0.22       345\n",
      "           4       0.70      0.74      0.72       544\n",
      "           5       0.44      0.09      0.15       227\n",
      "\n",
      "    accuracy                           0.54      5737\n",
      "   macro avg       0.50      0.40      0.41      5737\n",
      "weighted avg       0.52      0.54      0.51      5737\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sklearn.metrics.classification_report(y_test, rf_pred_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is seen that the F1 is the highest for the trips of type 0 and 2, and the lowest for type 3 and 5."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# n_samples is set as 5000 by default. \n",
    "# Setting n_samples as -1 means entire validation set\n",
    "# the score metric used is accuracy, aka number of records that are correctly predicted\n",
    "imp = rfpimp.importances(rf_clf, X_test, y_test, n_samples=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            Importance\n",
      "Feature               \n",
      "p_nearby_4    0.046191\n",
      "p_nearby_2    0.032595\n",
      "midnight      0.032595\n",
      "p_nearby_0    0.025972\n",
      "weekday       0.016036\n",
      "distance_m    0.008192\n",
      "mode          0.008018\n",
      "midday        0.007321\n",
      "evening       0.005578\n",
      "p_nearby_5    0.004881\n",
      "p_nearby_3    0.002789\n",
      "duration      0.001917\n",
      "p_nearby_1    0.001394\n",
      "temp          0.001046\n",
      "morning       0.000872\n",
      "afternoon     0.000872\n",
      "precip       -0.006624\n"
     ]
    }
   ],
   "source": [
    "print(imp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEHCAYAAAB1DlnCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZhcVbnv8W9BIhCCIAoKKiYCCUaQYBhluATXywUPgkwGZApwQUBwREQRjeI5h6s4oCicgExiIIJGIiLwLjAyyZAwgxCH5FwmQRQZTgIC6fvHWk0q1TV1eu/u6tTv8zw8qd619qq10+R91t613vVWenp6EBEpw0pDPQARWXEpwIhIaRRgRKQ0CjAiUhoFGBEpzYihHkAnq1QqqwNjgFeGeCgig2kksLCnp+d/BtqRAkxzY6655poHxo4dO9TjEBk0CxYsYLfddtsUeHCgfSnANPfK2LFjGTduXMuGF18yg0WLXxqEIQ2NUautyqEHf2yohyGDp5BZuwJMQRYtfold9jhgqIdRmhuuumyohyDDkB7yikhpFGBEpDTDIsCYhTFm4YEC+5tdZH8iUl/HBxizsHLB/e0DvFhknyJSX8uHvGZhDHANcDuwBTAfONQ9LqrTdiFwEfBh0nfp+7vHh83C6sAPgM3yZ05zj1fmvn8CrJ67ON493moWdga+CjwJTAQ+BIwwCxdVjwHYLp+zd/58A451j/s0uJbRwGeBo4GfNbrmSqVydG6zysyZMzn11FNb/TWJSB3tzmDGA9Pd4/uA54HjmrR9xj2+HzgbODEfOwW4wT1uBUwGvpWDztOA5fZTgO9X9bM1cIp7nNBkDDcA7zEL6+Q2hwMXNBnbacC3gT7BsVpPT8/0np6eLYF9p0yZ0qypiDTRboB51D3ekl9fAuzQpO0v8p/zSKtgAXYFTjYL9wBzgFWBDUiznHPNwv3A5cCEqn7ucI8Lmo3BPfaQZkAHm4W1SDOa39QblFmYCGzkHme1uFYRKUi762Bqd6VqtkvVy/nP16r6rwD7usdHqhuahWnAU8DmpGBXvVKtdplyozFcAPwqn3u5e3y1wbi2Aybl27gRwLpmYY573LnJtYjIALQ7g9nALGyXXx8I3NzPz7kWOMEsVADMwhb5+JrAk+5xCXAI0OyBbt0xuMcngCeALwMXNjrZPZ7tHtd3j2NIM7D5Ci4i5Wo3wPwBOMws3AesTXq+0h+nkW6H7stfD5+Wj/8o93sbMI6+s5Z2x/BT0i3UQ/0cl4iUqNJqT978Tc9V7nHTQRnRcjALZwF3u8cfF9lvpVIZ98gjjzyiXCTlInWT+fPnM378+PE9PT3zB9rXsM9FMgvzSDOfzw31WDqVgoMMlZYBxj0uBJaZvZiFWUDtHgZfcI/XFje09rjHSbXHzMLtwCo1hw9xj/eXNY5OTnZUoqIMleWawfQubOtU7nGboR6DiAyDVAERGb6GxTOYIh40m4VRpMV8G5LW6PzKPZ5czAhFpJ6On8EUnOx4hnvchJTPtL1Z2L3AvkWkRtckO+bx/ja//pdZuAt4R6vrF5Hl123JjgDkvKUPA9fXe79SqRxdqVTmAj+fOXNmq+5EpIGuSXbsZRZGAJcC33ePf6nXRtnUIsXopmTHXtOBP7rH77VoJyID1DXJjvlzv5E/89P9HL+ILIeuSXY0C+8gPQuaANxlFu4xC/+nn9chIv3Q7i3SEvd4TKtGeSuE3tdzgZ3z68XAx+u0/yPwvqpDX8zH55Ce1fS2W8iyz2dq7QCc22Jsj5Fu1URkkAyLhXbNdEqy46jVVu3YnJ9Rq6061EOQLqVkRxEpjZIdC6JsapG+Oj5VQESGryELMGZhT7NQN9nQLLQsjGYWzjMLzR78YhYuNAv71Tk+xixoByaRkg3ZQ173OBuYPYDzB/IV8xjgY8CMAfQhIi2UEmCqEiRvBrYF7iWtuP0asC5wEOlr5y3d4/FmYSzpH/uIfF5vPzsD04BnSA+a5wEHu8ceszAHONE9zjULRwJfIC24+yPwsns8Pnezk1n4LPA24CT3eAVwOimH6R7gIvf43TL+HkS6XZm3SBsBZ5LWuWxCmjHsQEqA/FJN2zOBs3My5F9r3tuCtPJ2AvBuYPvqN83C+sCppEBm+bOqrZc/dw9SYAE4GbjJPU6sF1yU7ChSjDIDzAL3eH9OA3gQuD4nJ97P0iTIXtuTEhAhJS9Wu8M9Ppb7uafOuVsDv3OP/3CPr5CSJqv90j0uyat839rOwJXsKFKMMgPMy1Wvl1T9vIT6t2aNEiir+6lOoOzVanVu9flaySsyiDplJe8twAGkbRgO6ue5dwDfNQtvAl4A9iXNkpp5AVijv4MUkf7plHUwnwI+YRbuJGU7t809Pg78B2nHvQg8BDzX4rT7gFfNwr1m4TPLMV4RaUPLyo7DgVkY7R5fzJtJzQLOd4+zBtpvfyo7nnPu+R29kveYo44Y6mHIMKHKjn1NMwuBtFPedcAvB3sASnYU6WuFCDDu8cTWrcql0qwifa0QAaYTXHzJDBYtfql1wxKo9rR0KgWYggxlNnWn3pqJdMq3SCKyAhoWM5giSsfmfiaRNgZfDbga+FReXSwiJej4GUzBpWPPBo4GNs7/7VZg3yJSo2tKx5qF9YA3usff558vBj5Ci0JtIrL8uql07NuBx6p+fiwf60PZ1CLF6KbSsfUSHes+f1E2tUgxuql07GPAO6p+fgdpgyoRKUnXlI51j08CL5iFbfM4DgWu7Od1iEg/dE3p2OxY4DzgT8Cf0QNekVJ1TenYqjENaC2NiLRvWCy0a0alY5UtLZ1LpWM7kJIXZUWh0rEFKTLZUcmLsqLo+FQBERm+hkWAMQtzzMKWLdpMNQtnDdaYRKS1YRFgRGR4Kqt07EnAS+7x+2bhu8Dm7nEXs/BBUr7QxaQysquQ1qMcnjftngR8BxhNKhc7NS+Q6+13JdLK3Ufd45fNwuGkr7afJCVAvpzbfZi08O4NwN9JpVD+BjwCfMA9/i33NR/Y1j0+U8bfg0i3K2sGcyOwY369JTDaLIwkrVe5n/SPP+Qkx7nAZ/P7PwD2y98MnQ/8e1WfI0gL6ubn4LIeKUhtTyoZW71O5mZS4NgCuIxUk3oJy9ZdCsC99YKLkh1FilHWOph5wCSzsAZpVnEXKdDsCMwmBYNbzAKkWcbvSdnSmwKej69Mmpn0+i/gZ+6xN+hsA8xxj38DMAszSauBIeUZzcxB6A1Ab9Lk+aT0gO8BR9Ag87qnp2c6ML1SqYybMmXKI/XaiEhrpQQY9/hK3hvmcOBWUqGzycCGpH/s7h4PrD7HLGwGPOget6O+W4HJZuHb7rE3KbJR0uUPgO+4x9l5b5lpeVyPmoWnzMIupADV3yqSItIPZT7kvZG0H8yNwE3AMaTi9bcB25uFjQDMwiizMI70fGSd3oRGszDSLLy3qr8fk7a5vDwXWLsd2NksvDnfXu1f1XZN4PH8+rCacZ1HulX6mXt8rbCrFZE+ygwwNwHrAb93j0+RtlO4Kd/STAUuzYmLtwGbuMd/AfsB/9cs3EsKRh+o7tA9fod0u/UT0jYP00i3VzEf7zWNFIhuIj0srjab9BC50cZUIlKQFaJ0bH/k9TTfdY87tmo7VKVjVepVhpJKxy4ns3AyacsGPXsRGQRdFWDc4+nA6WX0XWQ2tbKjZUXRVQGmEyhTWrqJAkxB2s2mVqa0dBPlIolIafo9g8mVAF4E3gjc6B5jg3YfIS3rb7VProisoJb7Fsk9fqVFk48AVwEKMCJdqq0AYxZOIZX5eJSUlTzPLFxIKkh/hVk4HdgTeBW4jlR8bU/gf5mFLwP7AruQ6kK/gbSr/yHucVHu53lSrtLbSImJV+TPPYlUzmQJ8Bv3eLJZ2BD4IbAOsAg4yj0+3GDcFwKLgU2Ad5FSFw4jFWi73T1ObetvSUSWS8tnMHkLhQNINaH3AbaqeX9tYG/gvbms6zfc462kFbOfd48T3eOfgV+4x63c4+akEiRHVnWzHinTeg/y18hmYXfSLGibfM43c9vpwAk54/pEUumTZt5ECm6fIRVo+y7wXmAzszCx3gnKphYpRjszmB2BWb3F7s3C7Jr3nyelAZxnFn5Nui2qZ1Oz8A1gLdJS/eoNwn+Zt1N4yCy8NR8LwAW9n+se/2EWRpPSBy7PGdfQd3PvWr9yjz25PO1TvRt/m4UHSaVt76k9QdnUIsVY3tKxr3OPr5qFrYEPkmY6x5NmDLUuBD7iHu81C1PJNZOyl6teV6r+rP3clYB/use6M48GevteUvM5S9DX9CKlaudr6huBvc3Canl/lw9Xv5lnFWu6x6uBTwO9//hfANaoaroG8GTOfG5nqf51wBFmYVT+nLXd4/PAArOwfz5WMQubt9GXiAyBlgHGPd4FzCTdSvyclCVdbQ3gqpwZ/TvSsw5IO8l93izcnR/MnkraYsGBug9laz73GtJznLlm4R7S8xZIwenInHH9ILBXq75EZGh0XTZ1f5SRTa1Mael0yqbuQO0mOyqRUbrJChFg8jqd/WsOX161f6+IDIEVIsDkQDKkwUTJjiJ9KdlRREqzwgcYs7DQLLxlqMch0o1W+AAjIkOnI5/BmIUxwDXkCo3AvaQqAF8D1iWthfkTqZDau0lJj0e7x/vMwpuBS0nJkHewdGUwZuFg4JOkhMvbgeNUukSkPJ08g9kIOBN4Hykb+mOkhMgTgS+Rgs3dOcHyS6R61wBfBW7OZWNnAxsAmIX3AFOA7XOqwWto82+RUnVygFngHu/PSZAPAte7xx5SbesxpGDzEwD3eAPwZrOwJrATqbAa7vHXwLO5vw8Ck4A788rgD5JmP30om1qkGB15i5TVJiZWJy2OIO09U6un5s9qFeAi9/jFVh+sbGqRYnTyDKaVG8m3OLn+9DM5GbL6+O6k/WAArgf2Mwvr5vfWNgvvGuxBi3ST4RxgpgFb5iTL01lag/prwE5m4S5gV+D/AeS9gb8MXJfPcdJGVyJSEiU7NqFkR+lGRSY7DucZjIh0uE5+yDusKJtapC/NYAaRysZKt9EMpiDtZFMrk1q6jWYwIlKajgwwZmFPs3Byg/debHD8QrOwX7kjE5H+6MhbJPc4m5RHJCLD2KAHmDYzpScAW7rH483CWGBGHus1Vf1UgB+QajAtYNms6a+QyqusBtwKfJyUd3S5e3x/brMxcFmuECkiJRiqW6RWmdLVzgTOdo9bAX+tOr43MB7YDDiKVPGx11m5TO2mpCCzRy5f+1xVudjDScXg+lCyo0gxhirAtMqUrrY9aX8XyNnT2U7Ape7xNff4BHBD1XuTzcLtuVzsLqRa1ADnAYebhZVJWzfMqDe4np6e6T09PVsC+06ZMmW5L1Kk2w1VgGmVKV2rUT5Dn+NmYVXgR8B+7nEz4Fygd3Xbz4HdgT2Aee7x7/0fuoi0qyO/RapxC6nmNSy7QdSNwAFmYWWzsB4wOR/vDSbP5LK2r3+z5B5fAq4FziY99xGREg2HAPMp4BNm4U5gzarjs4A/km6rziaVrcU9/pM0a7kf+CVwZ01/PyXNfK4rd9gi0nXZ1GbhRGBN93hqq7ZFZ1Mrk1qGA5WOXU5mYRawIenBb6HaSXZUoqN0m64KMO5x77L6VhKjSF9dFWDKdPElM1i0+KW67ymLWrqVAkxBmmVTK4tautVw+BapJbNwjFk4dKjHISLLWiFmMO7xnKEeg4j0NWgBpk7Z1vuAd7nHk/L7U4FJ7vGERiVe81YNZ5JW4i4G9nKPT5mFacCL7vEMszAnnzMZWAs40j3eZBZGkXKPNgH+QEpJ+IR7nFv+1Yt0p0G5RWpQtvVFYJ+qZlOAmS1KvK4O3OYeNyet5D2qwUeOcI9bA58mlZIFOA54NpeaPY1U5VFESjRYz2DqlW0dC/zFLGybC9aPJ6UFNCvx+i/gqvx6Hn0TI3v9ok6bHYDLANzjA6QZVF3KphYpxmDdItUt22oWjgQ+CjwMzHKPPXmfl0YlXl/JWdeQZjaNxv9ynTaVBm37UOlYkWIM1gymUdnWXwAfAQ4EZrZoO1A3k4IZZmECaR8ZESnRoASYRmVb3eOzwEOkh713NGtbwDB+BKyT+/wC6RbpuQL6FZEGuibZMW8yNdI9vmQWNiTNlMa5x381OqeoZEclOcpwomTH5TMK+K1ZGEl6HnNss+AiIgPXNQHGPb4AbFlW/82yqZVFLd1qhUgV6GRKdJRu1jUzmLI1SnZUoqN0M81gRKQ0CjAiUpphcYuUq0FelQupDaSfOaQ1NYvzoV3d49MDG52INNLxASavXynSQcqgFhkcLQNMVS3p24EtgPnAoe5xUZ22C4GLSHWhRwL7u8eHzcLqpDrSm+XPnOYer8x9/4SUJQ1wvHu81SzsTMqCfhKYCHwIGGEWLqoeA7BdPmfv/PlGWt9SnaXdb5VK5WjgaGCVmTNncuqpLQsQiEgd7T6DGQ9Mz1sdPE/a+qCRZ3KB+bNJtaYBTgFuyPWlJwPfykHnacBy+ynA96v62Ro4xT1OaDKGG4D3mIV1cpvDaV1Q7QKzcI9ZODUnVvah0rEixWg3wDzqHm/Jry8hbX3QSL2tEnYFTs7bL8whVV/cgDTLOTfXkL4cmFDVzx3ucUGzMeTM6p8AB5uFtUgzmt80GdtBuZzsjvm/Q5q0FZEBavcZTG3CUrMEpkZbJezrHpfZ+iDvRPcUsDkp2FVvy/8/bY7hAuBX+dzL3eOrjQbmHh/Pf75gFmaQZkkXN7kWERmAdmcwG5iF7fLrA0lbH/THtcAJvbckZmGLfHxN4En3uIQ0m2j2QLfuGNzjE8ATpAzsCxudbBZGmIW35NcjSdtuPtDP6xCRfmg3wPwBOCxvdbA26flKf5xGuh26zyw8kH+GtIXCYWbhNmAcfWct7Y7hp6RbqIeanL8KcG0+/x7gcVINaxEpScvtGopag1Ims3AWcLd7/HGR/fZnu4ZGhdeUiyTDjbZrqGIW5pFmPp8bynEoiIj01TLAuMeFwDKzl1xEfmxN0y+4x2uLG1p73GOf6gBm4XbSLVG1Q9zj/WWNo3oGo1mLSLJcM5gyi8gXwT1uM9ifWZ1NrQxqkUTJjiJSmmHxDKbAZMdrSMmOI4CbSJUdXxv4CEWkno6fwRSc7PjRXBVyU2AdYP8C+xaRGl2V7Ogen6+67jfQfEWyiAxQ1yU7moVr8+e+AFxRr41Kx4oUo9uSHXGP/5v0HGYVYJd6bZRNLVKMrkp27JWLr80G9iJVjhSREnRTsuNos7Befj2C9Fzn4X5eh4j0QzclO64OzM7n30t6DnNOP69DRPqh3VukJe7xmFaN3OOYqtdzgZ3z68XAx+u0/yPwvqpDX8zH55Ce1fS2W8iyz2dq7UCLzGj3+BSwVbM2IlKsYbHQrplOSXasLh2rUrEiiZIdC6ZER5GllOxYkN5kRyU6iizV8akCIjJ8DWqAMQvTzMKJrVu27Gcts3Bc1c/rm4W6q3JFZOh07Awmr1VpZC2q0hXc4xPucb/yRyUi/VH6t0hm4RRSYuKjwN+AeblG9InucW7e6X+uexxjFqYC/0ZKJVjdLOwJXAm8ibSO5svu8UrgdGDDnHrgwA/J2zmYhVVJa2S2BF4FPusef5v73hMYBWwIzHKPJ5V9/SLdrNQZjFmYBBxAyoDeh/bWoWwHHOYedyEt/987J0NOBr6dVwOfDPzZPU50j5+vOf8TALnA2oHARTnoQMrMnkLK6p5iFt5ZbwBKdhQpRtm3SDuSZgqL8lYJs9s4x93jP/LrCvAfefVtBN4OvLXF+TuQEiBxjw8D/01aJQxwvXt8zj2+BDwEvKteB0p2FCnGYDyDqZcY+WrVZ9euSqtOFziItDHUJPc4kZQY2WoVW91609nLVa+rkzFFpARlB5gbgb3NwmpmYQ3SRlQAC4HeBXLNHs6uCTztHl8xC5NZOuN4AVijyWceBGAWxpG2hXikQVsRKVGpAcY93gXMJFVS/DlpH1yAM4BjzcKtwFuadPFTYEuzMJcUNB7O/f4duMUsPGAWvlVzzo+AlfMeMzOBqe7xZURk0LWs7NjN+lPZ8Zxzz399Je8xRx0xCKMTKUeRlR07dh2MiAx/CjAF6c2mVia1yFIKMCJSGgWYgvRmU/fWpxYRBRgRKdGwWGhWYOnYfyflRb3JPY4uYmwi0ljHz2AKLh37K1JBNxEZBN1WOva23K7VZYtIAbqudGw7lE0tUoyuKx3bDmVTixSjK0vHisjg6JrSsSIy+LqpdCxm4Ztm4TFglFl4LM+gRKQkXVM6NvdzEqB9eEUGybBYaNdMp5WOVbKjyFIqHVsglY0VWZZKxxZESY4ifXV8qoCIDF8dHWBqS8SKyPDS0QGGmhKxIjK8dPq3SLUlYp8GPkp6gDvLPX61KhnzZmBb4F7S6t6vAesCB7nHO/Kalw1JxdveCXzTPbb8altEll+nz2BeLxFLCjAbk5IgJwKTzMJOud1GwJmkNTWbAB8jrY05EfhSVX/vI9W+3g74illYfzAuQqRbdXqAqbZr/u9u4C5SINk4v7fAPd6fUw4eJJWI7QHuZ2nCJcCV7nGxe3wG+C0N9oZRNrVIMTr9FqlaBfhP9/hf1QfzLVJ1YbUlVT8vYdlrbCtps6enZzowvVKpjJsyZYqqQoosp06fwVSXiL0WOMIsjAYwC283C+v2s7+9zMKqZuHNpDSGOwsbqYj00dEBprpELGDADOD3ef+YK2hcn7qRO4BfA7cBp+VMbBEpScffIrnH2rX3Z9Zp9noqg3ucWvV6IcumOcx3j0cXOT4RaayjZzAiMrx1/AymKO5xWpn9K4tapC/NYESkNF0zgymbsqlF+tIMRkRKs8IEGLNwdS5dIiIdYtjcIpmFEc1KkrjHDw3meESktdIDTDvZzsCfgPOBdwOLgKPd4305A3p9Uj7RM2bhOmBPYBQpM3pW3si7t2ztlsBoUvG1m4EPAI8De7nHxWZhK+DHpD18bwZ2d4/LbAcqIsUZrFukVtnOXwPuzmVhvwRcXHXuJFKA6F1wN5FUZnYzYIpZeGedz9sY+KF7fC/wT2DffPwC4Bj3uB2pMFxdSnYUKcZgBZhW2c47kErA4h5vAN5sFtbM587OZU96Xe8en3OPLwEPAe9q8Hn35NfzgDH5+cwa7vHWfHxGo8GqdKxIMQYrwLTKdq7UOac307m2GFt1X9XlaVu1qfcZIlKiTvkW6UbSsxjMws7AM+7x+SI/wD0+C7xgFrbNhw4osn8R6atTvkWaBlyQy8IuAg4r6XOOBM41C/9Dqhz5XEmfIyJApaen7p5LKySzMNo9vphfnwys5x4/1ah9pVIZ98gjjzwybty4ln2fc+75ABxz1BEFjVZkaMyfP5/x48eP7+npmT/QvjplBjNY/s0sfJF03f8NTC2qYyU7ivTVVTOY/urPDEZkRVHkDKZTHvIOexdfMoOLL2n4zbdIV+q2W6TSKJtapK9SA4xZ2B/4OvBX4DPA+u7x6jI/U0Q6R9m3SEcCx7nHyaQl/v1KSDQLmmGJDGOF/QM2C78klWRdlZR39DZSCsBYs3A1KR9oNbOwA/CfwFXAD0g5RSOAae7xSrMwlVR9cVVgdbPwddI6mWdIG3jPAw52jz1m4YPAGfn8O4Fj3ePLTY4vBC4CPgyMBPZ3jw8X9XcgIssqcgZzhHucRMpo/iTwQ2AuqTb0Z4CvADPd40T3OBM4BbjBPW4FTAa+ZRZWz31tBxzmHnfJP28BfBqYQMq43t4srApcCExxj71B6thGx6vG+Yx7fD9wNinZUkRKUmSA+aRZuJdUc+idLC3r2siuwMm5sP0c0oxlg/yeu8d/VLW9wz0+lpMl7yElSI4nJTX2fpV2EbBTk+O9fpH/nMeyZWVfp2xqkWIUcouU84cCsJ17XGQW5pACRjMVYF/3uExpVrOwDe0lODZKXmyV1NjbV6NESZWOFSlIUTOYNYFnc3DZhLSxVK3qMrCQSsGeYBYqAGZhi35+5sOkbRg2yj8fAvyuyXERGWRFBZhrgBE5WfE00m1Srd8CE8zCPWZhSm43Ergvl4Y9rT8fmPeDORy4PJeSXQKc0+j4cl6XiAyAUgWaULKjdCOlCojIsKCFbAVRNrVIXwowbbr4khlN841GrbYqhx78sYbvi3QjBZg2LVr8Ervs0XiXzRuuumwQRyMyPOgZjIiUZoUPMGbh1tatRKQMwyrAmIWV+3uOe/xAGWMRkdY65hlMVYnZ20nJjfOBQ0nF1c4n5S6dZRbuJCVSrkOqQHCUe3zYLLyVtKDu3bnLY93jrWbhRfc4OqczfB34Oylf6UbSVhJLBukSRbpOp81gxgPTcwnZ54Hj8vGX3OMO7vEyYDpwQs7cPhH4UW7zfeB37nFz4P2kCpK1tgY+R9oiYkNgn3qDULKjSDE6ZgaTPeoeb8mvLyFt+wAwE1LZEVJB+8vNQu85q+Q/dyHNeHCPr1G/5tEd7vEvua9LSfvVXFHbSMmOIsXotABTm7dQWz52JeCf7nFiwf2LSAk67RZpA7OwXX59IHBz9Zu5nOyCvNcvZqFiFjbPb19P3ljKLKxsFt5Yp/+tzcJYs7ASMKW2fxEpVqcFmD8Ah+Ws7LVJu87VOgg4Mm9u9SCwVz7+KWByzqCeB7y3zrm/B04HHgAWALOKHb6IVOu0W6Ql7vGYmmNjqn9wjwuA3WpPdI9PsTTYVB8fXfXjIvc4pYBxikgbOm0GIyIrkI6ZwbjHhaSqAWX1P4e0929/jFywYAEAz//zWS798fcaNlxtlVWYP3/A22eIDLn8//zIIvrShlNNVCqV1Um3aK+00fyjwM9KHdDQWtGvT5YaCUzu6ek5a6AdKcAUpFKpzO3p6dlyqMdRlhX9+mRZRf2+9QxGREqjACMipVGAKc70oR5AyVb065NlFfL71jMYESmNZjAiUpqOWQfTyczCbsCZwMrAee7x9Jr3K/n9D5H2qJnqHu9q59yhNsBrW0iq2Pka8Kp71LdMw0Abv/NNgAtI256c4h7PqHpvIf34nWsG00LeRe+HwO7ABOBAszChptnuwMb5v6PJOVRtnjtkBnJtVSa7x4kKLsNDm7/zf5C2SjmD+tr+nSvAtLY18Cf3+Bf3+C/gMvrmPO0FXOwee9zjbcBaZo3GPkUAAAFgSURBVGG9Ns8dSgO5NhmeWv7O3ePT7vFO2ltg2pQCTGtvBx6t+vmxfKydNu2cO5QGcm2Q9tO5zizMMwtHlzZKKdJA/5/s1+9cAaa1Sp1jtV+9NWrTzrlDaSDXBrC9e3w/abr9CbOwU5GDk1IM9P/Jfv3OFWBaewx4Z9XP7wCeaLNNO+cOpYFcG+6x98+nSXvrbF3aSKUoA/p/sr+/c32L1NqdwMZmYSzwOHAAUFsjdjZwvFm4DNgGeM49PmkW/tbGuUNpINe2OrCSe3whv96VVLVBOls7v/O6lud3rgDTgnt81SwcD1xL+lrvfPf4oFk4Jr9/DnA16WvcP5G+yj282blDcBl1DeTagLcCs/Lm6yOAGe7xmkG+BOmndn7nZuFtwFzgjcASs/Bp0jdOb6Gfv3Ot5BWR0ugZjIiURgFGREqjACMipVGAEZHSKMCISGkUYESkNAowIlIaBRgRKc3/B2Y8n1I6NlH3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x270.72 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "viz = rfpimp.plot_importances(imp)\n",
    "viz.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be seen that the most important variables are:\n",
    "- purpose_nearby_2 # shops\n",
    "- purpose_nearby_4 # education\n",
    "- mode\n",
    "- temp\n",
    "- purpose_nearby_0 # leisure\n",
    "- weekday\n",
    "- purpose_nearby_1 # food/drink\n",
    "- distance_m\n",
    "- purpose_nearby_3 # pick-up/drop-off\n",
    "- purpose_nearby_5 # health\n",
    "- midnight\n",
    "- duration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The hyperparameter of **n_estimators** will be tuned using *GridSearchCV*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/model_selection/_split.py:1978: FutureWarning: The default value of cv will change from 3 to 5 in version 0.22. Specify it explicitly to silence this warning.\n",
      "  warnings.warn(CV_WARNING, FutureWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters set found on development set:\n",
      "\n",
      "{'n_estimators': 200}\n",
      "\n",
      "Grid scores on development set:\n",
      "\n",
      "0.524 (+/-0.009) for {'n_estimators': 10}\n",
      "0.557 (+/-0.009) for {'n_estimators': 50}\n",
      "0.563 (+/-0.002) for {'n_estimators': 100}\n",
      "0.566 (+/-0.003) for {'n_estimators': 200}\n",
      "\n",
      "Detailed classification report:\n",
      "\n",
      "The model is trained on the full evaluation set.\n",
      "The scores are computed on the full evaluation set.\n",
      "\n",
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.54      0.68      0.60      2012\n",
      "           1       0.38      0.12      0.18       678\n",
      "           2       0.54      0.64      0.59      1931\n",
      "           3       0.45      0.15      0.23       345\n",
      "           4       0.71      0.74      0.72       544\n",
      "           5       0.36      0.08      0.13       227\n",
      "\n",
      "    accuracy                           0.55      5737\n",
      "   macro avg       0.50      0.40      0.41      5737\n",
      "weighted avg       0.53      0.55      0.52      5737\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# testing on different n_estimators\n",
    "list_n_estimators = [10, 50, 100, 200]\n",
    "\n",
    "# by default n_estimator in random forest is 100\n",
    "parameters = {'n_estimators': list_n_estimators}\n",
    "\n",
    "# use the default 5-fold cross validation. \n",
    "# About CV: For integer/None inputs, if the estimator is a classifier and y is either binary or multiclass, StratifiedKFold is used. \n",
    "# In all other cases, KFold is used.\n",
    "gscv_rfc = sklearn.model_selection.GridSearchCV(sklearn.ensemble.RandomForestClassifier(), parameters)\n",
    "gscv_rfc.fit(X_train, y_train)\n",
    "\n",
    "print(\"Best parameters set found on development set:\")\n",
    "print()\n",
    "print(gscv_rfc.best_params_)\n",
    "print()\n",
    "print(\"Grid scores on development set:\")\n",
    "print()\n",
    "means = gscv_rfc.cv_results_['mean_test_score']\n",
    "stds = gscv_rfc.cv_results_['std_test_score']\n",
    "for mean, std, params in zip(means, stds, gscv_rfc.cv_results_['params']):\n",
    "    print(\"%0.3f (+/-%0.03f) for %r\"\n",
    "          % (mean, std * 2, params))\n",
    "print()\n",
    "\n",
    "print(\"Detailed classification report:\")\n",
    "print()\n",
    "print(\"The model is trained on the full evaluation set.\")\n",
    "print(\"The scores are computed on the full evaluation set.\")\n",
    "print()\n",
    "y_true, y_pred = y_test, gscv_rfc.predict(X_test)\n",
    "print(sklearn.metrics.classification_report(y_true, y_pred))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/model_selection/_split.py:1978: FutureWarning: The default value of cv will change from 3 to 5 in version 0.22. Specify it explicitly to silence this warning.\n",
      "  warnings.warn(CV_WARNING, FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n",
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters set found on development set:\n",
      "\n",
      "{'max_depth': 10}\n",
      "\n",
      "Grid scores on development set:\n",
      "\n",
      "0.545 (+/-0.004) for {'max_depth': 10}\n",
      "0.528 (+/-0.010) for {'max_depth': 20}\n",
      "0.525 (+/-0.005) for {'max_depth': 50}\n",
      "0.522 (+/-0.011) for {'max_depth': 100}\n",
      "\n",
      "Detailed classification report:\n",
      "\n",
      "The model is trained on the full evaluation set.\n",
      "The scores are computed on the full evaluation set.\n",
      "\n",
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.54      0.68      0.60      2012\n",
      "           1       0.31      0.03      0.05       678\n",
      "           2       0.52      0.68      0.59      1931\n",
      "           3       0.42      0.09      0.14       345\n",
      "           4       0.71      0.70      0.70       544\n",
      "           5       0.42      0.07      0.12       227\n",
      "\n",
      "    accuracy                           0.54      5737\n",
      "   macro avg       0.49      0.37      0.37      5737\n",
      "weighted avg       0.51      0.54      0.50      5737\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.\n",
      "  \"10 in version 0.20 to 100 in 0.22.\", FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "# testing on different max_depth\n",
    "\n",
    "# by default n_estimator in random forest is 100\n",
    "parameters = {'max_depth': [10, 20, 50, 100]}\n",
    "\n",
    "# use the default 5-fold cross validation. \n",
    "# About CV: For integer/None inputs, if the estimator is a classifier and y is either binary or multiclass, StratifiedKFold is used. \n",
    "# In all other cases, KFold is used.\n",
    "gscv_rfc = sklearn.model_selection.GridSearchCV(sklearn.ensemble.RandomForestClassifier(), parameters)\n",
    "gscv_rfc.fit(X_train, y_train)\n",
    "\n",
    "print(\"Best parameters set found on development set:\")\n",
    "print()\n",
    "print(gscv_rfc.best_params_)\n",
    "print()\n",
    "print(\"Grid scores on development set:\")\n",
    "print()\n",
    "means = gscv_rfc.cv_results_['mean_test_score']\n",
    "stds = gscv_rfc.cv_results_['std_test_score']\n",
    "for mean, std, params in zip(means, stds, gscv_rfc.cv_results_['params']):\n",
    "    print(\"%0.3f (+/-%0.03f) for %r\"\n",
    "          % (mean, std * 2, params))\n",
    "print()\n",
    "\n",
    "print(\"Detailed classification report:\")\n",
    "print()\n",
    "print(\"The model is trained on the full evaluation set.\")\n",
    "print(\"The scores are computed on the full evaluation set.\")\n",
    "print()\n",
    "y_true, y_pred = y_test, gscv_rfc.predict(X_test)\n",
    "print(sklearn.metrics.classification_report(y_true, y_pred))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Applications/anaconda/envs/st-ds/lib/python3.7/site-packages/sklearn/model_selection/_split.py:1978: FutureWarning: The default value of cv will change from 3 to 5 in version 0.22. Specify it explicitly to silence this warning.\n",
      "  warnings.warn(CV_WARNING, FutureWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters set found on development set:\n",
      "\n",
      "{'max_depth': 100, 'n_estimators': 200}\n",
      "\n",
      "Grid scores on development set:\n",
      "\n",
      "0.544 (+/-0.010) for {'max_depth': 10, 'n_estimators': 10}\n",
      "0.554 (+/-0.006) for {'max_depth': 10, 'n_estimators': 50}\n",
      "0.555 (+/-0.008) for {'max_depth': 10, 'n_estimators': 100}\n",
      "0.555 (+/-0.011) for {'max_depth': 10, 'n_estimators': 200}\n",
      "0.531 (+/-0.014) for {'max_depth': 20, 'n_estimators': 10}\n",
      "0.556 (+/-0.006) for {'max_depth': 20, 'n_estimators': 50}\n",
      "0.564 (+/-0.002) for {'max_depth': 20, 'n_estimators': 100}\n",
      "0.564 (+/-0.002) for {'max_depth': 20, 'n_estimators': 200}\n",
      "0.522 (+/-0.010) for {'max_depth': 50, 'n_estimators': 10}\n",
      "0.556 (+/-0.006) for {'max_depth': 50, 'n_estimators': 50}\n",
      "0.562 (+/-0.004) for {'max_depth': 50, 'n_estimators': 100}\n",
      "0.561 (+/-0.000) for {'max_depth': 50, 'n_estimators': 200}\n",
      "0.527 (+/-0.007) for {'max_depth': 100, 'n_estimators': 10}\n",
      "0.557 (+/-0.006) for {'max_depth': 100, 'n_estimators': 50}\n",
      "0.562 (+/-0.002) for {'max_depth': 100, 'n_estimators': 100}\n",
      "0.565 (+/-0.002) for {'max_depth': 100, 'n_estimators': 200}\n",
      "\n",
      "Detailed classification report:\n",
      "\n",
      "The model is trained on the full evaluation set.\n",
      "The scores are computed on the full evaluation set.\n",
      "\n",
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.54      0.68      0.60      2012\n",
      "           1       0.38      0.11      0.18       678\n",
      "           2       0.54      0.64      0.59      1931\n",
      "           3       0.44      0.15      0.22       345\n",
      "           4       0.70      0.73      0.71       544\n",
      "           5       0.31      0.07      0.11       227\n",
      "\n",
      "    accuracy                           0.55      5737\n",
      "   macro avg       0.49      0.40      0.40      5737\n",
      "weighted avg       0.52      0.55      0.52      5737\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# testing n_estimators and max_depth\n",
    "\n",
    "# by default n_estimator in random forest is 100\n",
    "parameters = {'max_depth': [10, 20, 50, 100], 'n_estimators':[10, 50, 100, 200]}\n",
    "\n",
    "# use the default 5-fold cross validation. \n",
    "# About CV: For integer/None inputs, if the estimator is a classifier and y is either binary or multiclass, StratifiedKFold is used. \n",
    "# In all other cases, KFold is used.\n",
    "gscv_rfc = sklearn.model_selection.GridSearchCV(sklearn.ensemble.RandomForestClassifier(), parameters)\n",
    "gscv_rfc.fit(X_train, y_train)\n",
    "\n",
    "print(\"Best parameters set found on development set:\")\n",
    "print()\n",
    "print(gscv_rfc.best_params_)\n",
    "print()\n",
    "print(\"Grid scores on development set:\")\n",
    "print()\n",
    "means = gscv_rfc.cv_results_['mean_test_score']\n",
    "stds = gscv_rfc.cv_results_['std_test_score']\n",
    "for mean, std, params in zip(means, stds, gscv_rfc.cv_results_['params']):\n",
    "    print(\"%0.3f (+/-%0.03f) for %r\"\n",
    "          % (mean, std * 2, params))\n",
    "print()\n",
    "\n",
    "print(\"Detailed classification report:\")\n",
    "print()\n",
    "print(\"The model is trained on the full evaluation set.\")\n",
    "print(\"The scores are computed on the full evaluation set.\")\n",
    "print()\n",
    "y_true, y_pred = y_test, gscv_rfc.predict(X_test)\n",
    "print(sklearn.metrics.classification_report(y_true, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best hyperparameter found is **{'max_depth': 20, 'n_estimators': 200}**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Oversampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['mode', 'duration', 'distance_m', 'weekday', 'precip', 'temp',\n",
       "       'morning', 'midday', 'afternoon', 'evening', 'midnight', 'p_nearby_0',\n",
       "       'p_nearby_1', 'p_nearby_2', 'p_nearby_3', 'p_nearby_4', 'p_nearby_5'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define categorical variables\n",
    "list_categorical_var = ['mode', 'weekday', 'morning', 'midday', 'afternoon', 'evening', 'midnight']\n",
    "mark_cateogrical_var = [e in list_categorical_var for e in X_train.columns]\n",
    "rd_state = 42\n",
    "\n",
    "sm = SMOTENC(random_state=rd_state, categorical_features=mark_cateogrical_var)\n",
    "X_train_res, y_train_res = sm.fit_resample(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number per trip purpose in y_train\n",
      "0    6037\n",
      "2    5792\n",
      "1    2036\n",
      "4    1631\n",
      "3    1034\n",
      "5     681\n",
      "dtype: int64\n",
      "Number per trip purpose in y_train_res\n",
      "5    6037\n",
      "4    6037\n",
      "3    6037\n",
      "2    6037\n",
      "1    6037\n",
      "0    6037\n",
      "dtype: int64\n"
     ]
    }
   ],
   "source": [
    "# check number of trip purposes in y_train\n",
    "print(\"Number per trip purpose in y_train\")\n",
    "print(pd.Series(y_train).value_counts())\n",
    "print(\"Number per trip purpose in y_train_res\")\n",
    "print(pd.Series(y_train_res).value_counts())\n",
    "# check number of trip purposes in y_train_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy on training data:  0.9999418976236129\n",
      "Accuracy on testing data:  0.5478473069548545\n"
     ]
    }
   ],
   "source": [
    "rf, rf_accuracy_test, rf_pred_test = run_rf(X_train, X_test, y_train, y_test,\\\n",
    "                                                        n_estimators=200, max_depth=20)\n",
    "print(\"Accuracy on training data: \", rf_clf.score(X_train, y_train))\n",
    "print(\"Accuracy on testing data: \", rf_accuracy_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.54      0.67      0.60      2012\n",
      "           1       0.40      0.11      0.17       678\n",
      "           2       0.54      0.65      0.59      1931\n",
      "           3       0.44      0.14      0.21       345\n",
      "           4       0.71      0.74      0.72       544\n",
      "           5       0.38      0.08      0.13       227\n",
      "\n",
      "    accuracy                           0.55      5737\n",
      "   macro avg       0.50      0.40      0.40      5737\n",
      "weighted avg       0.53      0.55      0.51      5737\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sklearn.metrics.classification_report(y_test, rf_pred_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy on training data:  0.6590193804869968\n",
      "Accuracy on testing data:  0.5291964441345651\n"
     ]
    }
   ],
   "source": [
    "# compare the model using resampled data and \n",
    "rf_res, rf_res_accuracy_test, rf_res_pred_test = run_rf(X_train_res, X_test, y_train_res, y_test,\\\n",
    "                                                        n_estimators=200, max_depth=20)\n",
    "# print results\n",
    "print(\"Accuracy on training data: \", rf_clf.score(X_train_res, y_train_res))\n",
    "print(\"Accuracy on testing data: \", rf_res_accuracy_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The encoding of **purpose** and **mode** is as follows. Note that the trips of home and work have been removed.\n",
    "\n",
    "- purpose: {'leisure': 0, 'food_drink': 1, 'shops': 2, 'pick_up_drop_off': 3, 'education': 4, 'health': 5}\n",
    "- mode: {'walking': 0, 'public_transport': 1, 'car': 2, 'cycling': 3, 'public_transport, car': 4, 'other': 5}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.58      0.61      0.60      2012\n",
      "           1       0.32      0.29      0.30       678\n",
      "           2       0.58      0.54      0.56      1931\n",
      "           3       0.29      0.31      0.30       345\n",
      "           4       0.66      0.78      0.71       544\n",
      "           5       0.23      0.22      0.23       227\n",
      "\n",
      "    accuracy                           0.53      5737\n",
      "   macro avg       0.44      0.46      0.45      5737\n",
      "weighted avg       0.53      0.53      0.53      5737\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sklearn.metrics.classification_report(y_test, rf_res_pred_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The oversampling doesn't lead to increase on the prediction accuracy. More predictor variables are needed for this task."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Export model results\n",
    "#### Functions in this section:\n",
    "1. `model_setup` – Randomises the stratified train/test split and Trip IDs used  \n",
    "2. `merge_predictions_to_data` – Code to join predictions of testing set back to original data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Globals for columns to keep in the X data\n",
    "BASIC_VARIABLES = ['mode', 'duration', 'distance_m', 'weekday', 'precip',\n",
    "       'temp', 'morning', 'midday', 'afternoon', 'evening', 'midnight']\n",
    "TEMPORARY_VARIABLES = ['id_trip','endx','endy'] # needed for nearby purpose computation \n",
    "\n",
    "\n",
    "def model_setup(X, y, rd_state):\n",
    "    \"\"\"\n",
    "        Sets up the model by carrying out the train/test splits as well as calculating the nearby purposes columns\n",
    "        :param X: (dataframe) Independent variables\n",
    "        :param y: (dataframe) Dependent variable\n",
    "        :param rd_state: (int) Random state used for stratified train/test split\n",
    "    \"\"\"\n",
    "    ## extract all variables needed based on global column list (see globals)\n",
    "    X = X[BASIC_VARIABLES + TEMPORARY_VARIABLES]\n",
    "    \n",
    "    ## perform stratified train/test split\n",
    "    X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(\n",
    "        X, y, stratify=y, random_state=rd_state)\n",
    "    \n",
    "    ## compute nearby purposes\n",
    "    X_train, X_test = compute_nearby_purposes(y_train, X_train, X_test, 200)\n",
    "    \n",
    "    ## preserve Trip IDs of training and testing data\n",
    "    X_train_IDs = X_train['id_trip']\n",
    "    X_test_IDs = X_test['id_trip']\n",
    "    \n",
    "    # remove columns which are not needed\n",
    "    X_train = X_train.drop(columns=TEMPORARY_VARIABLES, axis=1)\n",
    "    X_test = X_test.drop(columns=TEMPORARY_VARIABLES, axis=1)\n",
    "    \n",
    "    # flatten y train data\n",
    "    y_train = np.ravel(y_train)\n",
    "    y_test = np.ravel(y_test)\n",
    "    \n",
    "    return X_train, X_test, y_train, y_test, X_train_IDs, X_test_IDs\n",
    "\n",
    "\n",
    "def merge_predictions_to_data(original_data, predictions, X_test, y_test, test_ids):\n",
    "    \"\"\"\n",
    "        Merge the predictions made by the model on the testing set back to the original data by merging on the trip ID\n",
    "        :param original_data: (dataframe) Original dataframe from which train-test split is carried out\n",
    "        :param predictions: (array) Predictions made on the testing set\n",
    "        :param test_ids: (array) trip IDs for each row in the testing set\n",
    "        :return original_data_with_preds: (dataframe) New dataframe containing only the testing set trip IDs and with predictions for these\n",
    "    \"\"\"\n",
    "    ## make predictions into dataframe\n",
    "    # apply series to the array of predictions\n",
    "    predictions = pd.Series(predictions)\n",
    "    # make a dataframe of ID, purpose and prediction\n",
    "    prediction_data = pd.DataFrame([test_ids, y_test, predictions]).T\n",
    "    prediction_data.columns = ['id_trip','purpose','prediction']\n",
    "    # join back to X_test\n",
    "    X_test = pd.concat([X_test, prediction_data],axis=1)\n",
    "    \n",
    "    ## merge to original\n",
    "    # get columns to perform merge on (to miss out duplicates)\n",
    "    cols_to_use = list(prediction_data.columns.difference(X_17.columns))\n",
    "    cols_to_use += ['id_trip', 'p_nearby_0','p_nearby_1','p_nearby_2','p_nearby_3','p_nearby_4','p_nearby_5']\n",
    "\n",
    "    # perform inner join\n",
    "    original_data_with_preds = original_data.merge(X_test[cols_to_use], on='id_trip', how='inner')\n",
    "    return original_data_with_preds\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "605\n"
     ]
    }
   ],
   "source": [
    "# set a random seed for train test split\n",
    "random.seed(datetime.datetime.now())\n",
    "rd_state = random.randint(0,1001)\n",
    "print(rd_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 36.4 s, sys: 634 ms, total: 37.1 s\n",
      "Wall time: 37.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# create new set of train/test data\n",
    "X_train, X_test, y_train, y_test, X_train_IDs, X_test_IDs = model_setup(X_17, y_17, rd_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy on testing data: 0.5492417639881472\n"
     ]
    }
   ],
   "source": [
    "# carry out modelling\n",
    "rf, rf_accuracy_test, rf_preds = run_rf(X_train, X_test, y_train, y_test,\\\n",
    "                                                        n_estimators=200, max_depth=20)\n",
    "print('accuracy on testing data:', rf_accuracy_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mode                    3\n",
      "duration             3603\n",
      "distance_m        2862.76\n",
      "weekday                 1\n",
      "precip        0.000750916\n",
      "temp              16.0571\n",
      "morning             False\n",
      "midday              False\n",
      "afternoon           False\n",
      "evening             False\n",
      "midnight             True\n",
      "p_nearby_0       0.705882\n",
      "p_nearby_1       0.117647\n",
      "p_nearby_2       0.176471\n",
      "p_nearby_3              0\n",
      "p_nearby_4              0\n",
      "p_nearby_5              0\n",
      "Name: 0, dtype: object\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id_trip</th>\n",
       "      <th>mode</th>\n",
       "      <th>duration</th>\n",
       "      <th>distance_m</th>\n",
       "      <th>weekday</th>\n",
       "      <th>precip</th>\n",
       "      <th>temp</th>\n",
       "      <th>morning</th>\n",
       "      <th>midday</th>\n",
       "      <th>afternoon</th>\n",
       "      <th>evening</th>\n",
       "      <th>midnight</th>\n",
       "      <th>startx</th>\n",
       "      <th>starty</th>\n",
       "      <th>endx</th>\n",
       "      <th>endy</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>17285</th>\n",
       "      <td>435174</td>\n",
       "      <td>3</td>\n",
       "      <td>3603</td>\n",
       "      <td>2862.764494</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000751</td>\n",
       "      <td>16.05711</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>7.624681e+06</td>\n",
       "      <td>1.246457e+06</td>\n",
       "      <td>7.626269e+06</td>\n",
       "      <td>1.247517e+06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       id_trip  mode  duration   distance_m  weekday    precip      temp  \\\n",
       "17285   435174     3      3603  2862.764494        1  0.000751  16.05711   \n",
       "\n",
       "       morning  midday  afternoon  evening  midnight        startx  \\\n",
       "17285        0       0          0        0         1  7.624681e+06   \n",
       "\n",
       "             starty          endx          endy  \n",
       "17285  1.246457e+06  7.626269e+06  1.247517e+06  "
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## DATA CHECK: to see if the ids still match up between the X_test and the original data\n",
    "print(X_test.iloc[0])\n",
    "X_all.loc[X_all.id_trip == X_test_IDs[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge predictions back to data\n",
    "data_with_preds = merge_predictions_to_data(X_all, rf_preds, X_test, y_test, X_test_IDs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id_trip</th>\n",
       "      <th>mode</th>\n",
       "      <th>duration</th>\n",
       "      <th>distance_m</th>\n",
       "      <th>weekday</th>\n",
       "      <th>precip</th>\n",
       "      <th>temp</th>\n",
       "      <th>morning</th>\n",
       "      <th>midday</th>\n",
       "      <th>afternoon</th>\n",
       "      <th>...</th>\n",
       "      <th>endx</th>\n",
       "      <th>endy</th>\n",
       "      <th>prediction</th>\n",
       "      <th>purpose</th>\n",
       "      <th>p_nearby_0</th>\n",
       "      <th>p_nearby_1</th>\n",
       "      <th>p_nearby_2</th>\n",
       "      <th>p_nearby_3</th>\n",
       "      <th>p_nearby_4</th>\n",
       "      <th>p_nearby_5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>150744</td>\n",
       "      <td>2</td>\n",
       "      <td>862</td>\n",
       "      <td>9935.922336</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>16.910884</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>7.631864e+06</td>\n",
       "      <td>1.250415e+06</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>0.766667</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>257884</td>\n",
       "      <td>3</td>\n",
       "      <td>110</td>\n",
       "      <td>347.261519</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000062</td>\n",
       "      <td>26.659408</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>7.630828e+06</td>\n",
       "      <td>1.250138e+06</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0.250000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.250000</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>410129</td>\n",
       "      <td>2</td>\n",
       "      <td>569</td>\n",
       "      <td>4605.407008</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000099</td>\n",
       "      <td>25.860051</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>7.658376e+06</td>\n",
       "      <td>1.256748e+06</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>45863</td>\n",
       "      <td>2</td>\n",
       "      <td>1037</td>\n",
       "      <td>16426.302328</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000099</td>\n",
       "      <td>22.094541</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>7.628635e+06</td>\n",
       "      <td>1.250061e+06</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.357143</td>\n",
       "      <td>0.071429</td>\n",
       "      <td>0.357143</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.142857</td>\n",
       "      <td>0.071429</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>452030</td>\n",
       "      <td>2</td>\n",
       "      <td>1008</td>\n",
       "      <td>6500.370675</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>18.627913</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>7.615643e+06</td>\n",
       "      <td>1.247518e+06</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 24 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   id_trip  mode  duration    distance_m  weekday    precip       temp  \\\n",
       "0   150744     2       862   9935.922336        1  0.000000  16.910884   \n",
       "1   257884     3       110    347.261519        1  0.000062  26.659408   \n",
       "2   410129     2       569   4605.407008        1  0.000099  25.860051   \n",
       "3    45863     2      1037  16426.302328        1  0.000099  22.094541   \n",
       "4   452030     2      1008   6500.370675        1  0.000000  18.627913   \n",
       "\n",
       "   morning  midday  afternoon  ...          endx          endy  prediction  \\\n",
       "0        1       0          0  ...  7.631864e+06  1.250415e+06           2   \n",
       "1        0       0          0  ...  7.630828e+06  1.250138e+06           2   \n",
       "2        0       0          0  ...  7.658376e+06  1.256748e+06           2   \n",
       "3        0       0          0  ...  7.628635e+06  1.250061e+06           0   \n",
       "4        0       1          0  ...  7.615643e+06  1.247518e+06           2   \n",
       "\n",
       "   purpose  p_nearby_0  p_nearby_1  p_nearby_2  p_nearby_3  p_nearby_4  \\\n",
       "0        3    0.133333    0.033333    0.766667    0.033333    0.033333   \n",
       "1        2    0.250000    0.000000    0.250000    0.500000    0.000000   \n",
       "2        3    0.000000    0.000000    0.000000    0.000000    0.000000   \n",
       "3        2    0.357143    0.071429    0.357143    0.000000    0.142857   \n",
       "4        5    0.000000    0.000000    1.000000    0.000000    0.000000   \n",
       "\n",
       "   p_nearby_5  \n",
       "0    0.000000  \n",
       "1    0.000000  \n",
       "2    0.000000  \n",
       "3    0.071429  \n",
       "4    0.000000  \n",
       "\n",
       "[5 rows x 24 columns]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# view data\n",
    "data_with_preds.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# append times for start and end of trip\n",
    "start_end_times = gpd.read_file('../../Data/mtl_trajet/mtl_trajet_2017_final.shp')[['id_trip','starttime','endtime']]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_with_preds = data_with_preds.merge(start_end_times, on='id_trip')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id_trip</th>\n",
       "      <th>mode</th>\n",
       "      <th>duration</th>\n",
       "      <th>distance_m</th>\n",
       "      <th>weekday</th>\n",
       "      <th>precip</th>\n",
       "      <th>temp</th>\n",
       "      <th>morning</th>\n",
       "      <th>midday</th>\n",
       "      <th>afternoon</th>\n",
       "      <th>...</th>\n",
       "      <th>prediction</th>\n",
       "      <th>purpose</th>\n",
       "      <th>p_nearby_0</th>\n",
       "      <th>p_nearby_1</th>\n",
       "      <th>p_nearby_2</th>\n",
       "      <th>p_nearby_3</th>\n",
       "      <th>p_nearby_4</th>\n",
       "      <th>p_nearby_5</th>\n",
       "      <th>starttime</th>\n",
       "      <th>endtime</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>150744</td>\n",
       "      <td>2</td>\n",
       "      <td>862</td>\n",
       "      <td>9935.922336</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>16.910884</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>0.133333</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>0.766667</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2017-09-18 10:02:50</td>\n",
       "      <td>2017-09-18 10:17:12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>257884</td>\n",
       "      <td>3</td>\n",
       "      <td>110</td>\n",
       "      <td>347.261519</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000062</td>\n",
       "      <td>26.659408</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0.250000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.250000</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2017-09-18 20:50:10</td>\n",
       "      <td>2017-09-18 20:52:00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>410129</td>\n",
       "      <td>2</td>\n",
       "      <td>569</td>\n",
       "      <td>4605.407008</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000099</td>\n",
       "      <td>25.860051</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2017-09-18 21:20:24</td>\n",
       "      <td>2017-09-18 21:29:53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>45863</td>\n",
       "      <td>2</td>\n",
       "      <td>1037</td>\n",
       "      <td>16426.302328</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000099</td>\n",
       "      <td>22.094541</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.357143</td>\n",
       "      <td>0.071429</td>\n",
       "      <td>0.357143</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.142857</td>\n",
       "      <td>0.071429</td>\n",
       "      <td>2017-09-19 00:40:14</td>\n",
       "      <td>2017-09-19 00:57:31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>452030</td>\n",
       "      <td>2</td>\n",
       "      <td>1008</td>\n",
       "      <td>6500.370675</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>18.627913</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2017-09-19 11:48:18</td>\n",
       "      <td>2017-09-19 12:05:06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 26 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   id_trip  mode  duration    distance_m  weekday    precip       temp  \\\n",
       "0   150744     2       862   9935.922336        1  0.000000  16.910884   \n",
       "1   257884     3       110    347.261519        1  0.000062  26.659408   \n",
       "2   410129     2       569   4605.407008        1  0.000099  25.860051   \n",
       "3    45863     2      1037  16426.302328        1  0.000099  22.094541   \n",
       "4   452030     2      1008   6500.370675        1  0.000000  18.627913   \n",
       "\n",
       "   morning  midday  afternoon  ...  prediction  purpose  p_nearby_0  \\\n",
       "0        1       0          0  ...           2        3    0.133333   \n",
       "1        0       0          0  ...           2        2    0.250000   \n",
       "2        0       0          0  ...           2        3    0.000000   \n",
       "3        0       0          0  ...           0        2    0.357143   \n",
       "4        0       1          0  ...           2        5    0.000000   \n",
       "\n",
       "   p_nearby_1  p_nearby_2  p_nearby_3  p_nearby_4  p_nearby_5  \\\n",
       "0    0.033333    0.766667    0.033333    0.033333    0.000000   \n",
       "1    0.000000    0.250000    0.500000    0.000000    0.000000   \n",
       "2    0.000000    0.000000    0.000000    0.000000    0.000000   \n",
       "3    0.071429    0.357143    0.000000    0.142857    0.071429   \n",
       "4    0.000000    1.000000    0.000000    0.000000    0.000000   \n",
       "\n",
       "             starttime              endtime  \n",
       "0  2017-09-18 10:02:50  2017-09-18 10:17:12  \n",
       "1  2017-09-18 20:50:10  2017-09-18 20:52:00  \n",
       "2  2017-09-18 21:20:24  2017-09-18 21:29:53  \n",
       "3  2017-09-19 00:40:14  2017-09-19 00:57:31  \n",
       "4  2017-09-19 11:48:18  2017-09-19 12:05:06  \n",
       "\n",
       "[5 rows x 26 columns]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_with_preds.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save data \n",
    "data_with_preds.to_csv('../../Data/model_outputs/predictions_17.csv',index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# POI_RF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One vs many\n",
    "see: https://stackoverflow.com/questions/56090541/how-to-plot-precision-and-recall-of-multiclass-classifier/56092736"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 37.5 s, sys: 634 ms, total: 38.1 s\n",
      "Wall time: 38.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# create new set of train/test data\n",
    "X_train, X_test, y_train, y_test, X_train_IDs, X_test_IDs = model_setup(X_17, y_17, rd_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "sklearn.preprocessing.label_binarize(y_train, classes=[0,1,2,3,4,5])\n",
    "y_train = sklearn.preprocessing.label_binarize(y_train, classes=[0,1,2,3,4,5])\n",
    "y_test = sklearn.preprocessing.label_binarize(y_test, classes=[0,1,2,3,4,5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn.multiclass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = sklearn.multiclass.OneVsRestClassifier(sklearn.ensemble.RandomForestClassifier(n_estimators=200,\n",
    "                             max_depth=20,\n",
    "                             random_state=42))\n",
    "clf.fit(X_train, y_train)\n",
    "\n",
    "y_score = clf.predict_proba(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3780721631514729"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydeVhUZfvHP88Mw74ooqKi4L4iuG9lLplbaWml2evW22L1VlZaVlZWr2m/t301bXXLMttM08o0NbeUVMRdREFEBWRfZph5fn+cYWBgwFEYAXk+18XFnGc79zkM5z7P9r2FlBKFQqFQ1F50VW2AQqFQKKoW5QgUCoWilqMcgUKhUNRylCNQKBSKWo5yBAqFQlHLUY5AoVAoajnKESiqFUKIu4UQvzpRboEQ4vmrYVNNQgghhRCtrJ+/EEL8t6ptUlR/3KraAIWiOFLKZcAyJ8pNuwrmKBS1AtUjUFQ6Qgj1glEO19L9uZaupTajHIHCKYQQcUKIZ4QQB4UQF4UQnwshPK15A4QQCUKIp4UQScDn1vSbhRB7hRBpQohtQojOxdprKoT4TghxQQiRIoR435o+RQix1fpZCCHeEkKcF0KkCyH2CyE6WfPshj2EEPcJIY4LIVKFED8JIRoXy5NCiGlCiGNW2z8QQggH19hYCJErhAgsltZFCJEshDAIIVoJIf602pIshPjayXs3RQjxl/VaUoE51vR7hBCHrDatF0KEFqvTUQjxm/V6zgkhnrWm9xRCbLfe07NCiPeFEO7O2OHArvus58+0/l27FrtfrYqVs91rR39raxs3FyvvZr0/he31tv7904QQ+4QQA67EXoXrUI5AcTncDQwFWgJtgNnF8oKBQCAUuN/6EPgMeACoB3wM/CSE8BBC6IGfgVNAGNAEWOHgfDcB/a3nqgOMA1JKFhJCDALmAXcCjaztlmzvZqAHEGEtN7RkO1LKRGA7MLZY8gTgWymlCXgF+BWoC4QA7zmwuSx6AbFAA2CuEOJW4FlgDFAf2AJ8Zb0eP+B3YB3QGGgFbLC2YwYeB4KAPsBg4KHLsAPrOe5Ac0iTAH9gFA7ubRnY/a2tdt9VLH8okCyljBJCNAHWAP+11pkBrBJC1L9cmxWuQzkCxeXwvpQyXkqZCszF/p/fArwopcyXUuYC9wEfSyl3SinNUsovgXygN9AT7QE3U0qZLaXMk1JudXA+E+AHtAOElPKQlPKsg3J3A59JKaOklPnAM0AfIURYsTLzpZRpUsrTwEYgsoxrXF54XdZew3hrWqE9oUDjcmwui0Qp5XtSygLr/XkAmGe9pgLgVSDS2iu4GUiSUr5hPU+mlHIngJRyj5Ryh7WdODQHe8Nl2FHIvcD/SSn/lhrHpZSnnKxb8m+9HBglhPC25k+g6J79C1grpVwrpbRIKX8DdgMjrsBmhYtQjkBxOcQX+3wK7WFeyAUpZV6x41DgSetwQJoQIg1oaq3TFDhlfQCWiZTyD+B94APgnBBioRDC30HRxlZ7Cutlob3dNilWJqnY5xzAt4zTfovmRBqj9UYk2ts6wFOAAHYJIWKEEPeUZ38J4kschwLvFLs3qda2m6DdnxOOGhFCtBFC/CyESBJCZKA5kKDLsKOQMs/hBHZ/aynlceAQcIvVGYyiyBGEAneU+B5ch9ZzU1QTlCNQXA5Ni31uBiQWOy4pYxsPzJVS1in24y2l/Mqa10w4MdEopXxXStkN6Ig2RDTTQbFEtAcOAEIIH7ThqDPOXFSJ86WhDf/cifZm+5W0SvRKKZOklPdJKRujvdF/WHws/VJNlziOBx4ocX+8pJTbrHkty2jnI+Aw0FpK6Y82vFRqvsMJyjtHDuBd7Di4RL4jyeLC4aHRwEGrcyg8z5IS1+kjpZx/BTYrXIRyBIrL4WEhRIh1MvVZoLzJ0kXANCFEL+ukr48QYqR1/HsXcBaYb033FEL0K9mAEKKHtb4ByAby0MbIS7IcmCqEiBRCeKC9Je+0Dp1cCcvRxs7HUvRmixDiDiFEiPXwItoD0ZE9zrAAeEYI0dHadoB13B60+ZNgIcR065yKnxCilzXPD8gAsoQQ7YAHr/D8nwAzhBDdrH+fVsUmq/cCE4QQeiHEMJwbelqBNqfzIMXuGbAUracw1Nqep3XCOcRhK4oqQTkCxeWwHO1tOdb6U+ZmJSnlbrR5gvfRHprHgSnWPDNwC9ok6GkgAW0iuCT+aA7lItrQTwrwuoNzbQCeB1ahOZiWaGP7V8pPQGvgnJRyX7H0HsBOIUSWtcxjUsqTANahorudPYGU8nvgNWCFdYjnADDcmpcJDEG7R0nAMWCgteoMtJ5KJtq9cWrlkoPzr0Sb51lubesHtMlcgMes505Dm3/5wYn2zqJNtPctbpOUMh6tl/AscAGthzAT9eypVggVmEbhDEKIOOBeKeXvVW2LQqGoXJRXVigUilqOcgQKhUJRy1FDQwqFQlHLUT0ChUKhqOXUOMGooKAgGRYWVtVmKBQKRY1iz549yVJKh9IeNc4RhIWFsXv37qo2Q6FQKGoUQogyJUTU0JBCoVDUcpQjUCgUilqOcgQKhUJRy6lxcwQKhUJRiMlkIiEhgby8vEsXriV4enoSEhKCwWBwuo5yBAqFosaSkJCAn58fYWFhiNJB52odUkpSUlJISEigefPmTtdz2dCQEOIzoYUYPFBGvhBCvCu08IL7C8PaKRQKhbPk5eVRr1495QSsCCGoV6/eZfeQXDlH8AUwrJz84WgKj63Rwt195EJbFArFNYpyAvZcyf1w2dCQlHJziVCBJRkNLLYG/dghhKgjhGhURijCCvPR1PvIz83HU98CT0NLPAx6ANzq10cfGHiJ2pePb11P+t3eCp8Aj0pvW6FQKCqTqpwjaIJ9+L4Ea1opRyCEuB+t10CzZs2u6GTmfAtmmUq+IRKzZxOyCzNSgdT0K2qzfNJJTczm1ie64Onj/KSNq5BSYi6wYDZZKDBZ8PQxoHdTi8YUClcwZ84cfH19mTFjRqW3vWfPHqZMmUJubi4jRozgnXfeqXCvqCodgSPLHSrgSSkXAgsBunfvfkUqeQHN/Tl//BwJvhvxlonc1qkvOV98jntIExrPn3clTZaJxSzZvOIoKWeyWPPBPkY91gWDh77cOlJKzCYLxjwzxrwCTHlmjLkFGPMKtLRin035ZsxGMwUFFsxGi/bbZKHAaMFcYKHAaNZ+m4oe/GaTxf5+NPBi/OyeuLmXb5dCoahePPjggyxcuJDevXszYsQI1q1bx/DhwyvUZlU6ggTsY+CGYB8Dt1Jp1KIN548fI65+MgbRnEZ92pL6znE8GkDj1nUr/Xy3PBrJd6/vISk2g18W7KfbsDByMo3kZhrJyTCSm2EkJ9Nkd1xQ4mFd2ejcBG4GPeYCC+nnc9m/MYGuQ0MvXVGhUJTJ4sWLef311xFC0LlzZ5YsWWKXv2jRIhYuXIjRaKRVq1YsWbIEb29vVq5cyUsvvYRerycgIIDNmzcTExPD1KlTMRqNWCwWVq1aRevWrW1tnT17loyMDPr06QPApEmT+OGHH2q0I/gJ+I8QYgXQC0h31fxAVeAX6Mnox7rw3et7iD90kfhDFy9ZR+cm8PByw+DphrunHndPN9y9in/WY/B0w+Chx82gw82gQ2/QW39rP2623/pin3Xo3XQIndYJiz+Yyk/v7mX3L3G06tYA/yAvV98OhcLlhM1a45J24+aPLDMvJiaGuXPn8tdffxEUFERqamqpMmPGjOG+++4DYPbs2Xz66ac88sgjvPzyy6xfv54mTZqQlpYGwIIFC3jssce4++67MRqNmM32IbHPnDlDSEhRuOeQkBDOnDlT4Wt0mSMQQnwFDACChBAJwIuAAUBKuQBYC4xAi2WbA0x1lS1VRZ2G3ox6LJK/vj2OucCCt587Xn7uePm74+3vbj022I4NHvqrsgKiaYdAmravS/yhiyyZvZ3wG5rQ/662Lj+vQnGt8ccff3D77bcTFBQEQKCDhScHDhxg9uzZpKWlkZWVxdChQwHo168fU6ZM4c4772TMmDEA9OnTh7lz55KQkMCYMWPsegOgDSGXpDKeGa5cNXTXJfIl8LCrzl9dCArxY/T0LlVtRin63d6aFa/sAuDIrnOkXcglN9NIXpYJdy83bnuya7WY5FYonKW8N3dXIaW85IN4ypQp/PDDD0RERPDFF1+wadMmQHv737lzJ2vWrCEyMpK9e/cyYcIEevXqxZo1axg6dCiffPIJgwYNsrUVEhJCQkKC7TghIYHGjRtX+DrUspFaSr0mvvzrld4AGHMLiD+YSnJ8FlkX80lNzGbp89v54c0ols/ZwSdPbubLZ/8iOy0fAGmRmtPINlXlJSgUVc7gwYP55ptvSElJAXA4NJSZmUmjRo0wmUwsW7bMln7ixAl69erFyy+/TFBQEPHx8cTGxtKiRQseffRRRo0axf79++3aatSoEX5+fuzYsQMpJYsXL2b06NEVvg4lMVGLCajvzS2PRJCTacTLVxum+nP5Ec6fyiQ/p4AzR9NsZfOzC/hi1l94B7iTm2lCWiRCwOjHu+BXzxO/QE+1sUdR6+jYsSPPPfccN9xwA3q9ni5duvDFF1/YlXnllVfo1asXoaGhhIeHk5mZCcDMmTM5duwYUkoGDx5MREQE8+fPZ+nSpRgMBoKDg3nhhRdKnfOjjz6yLR8dPnx4hSeKoQbGLO7evbu8ksA0v3/6Eft+XcOODikYRA/eGvUAqeNvx6N1K1qsXu0CS2smmal5xO1PxsPbTZu78HNn188nif3nQrn1etzcnJ43X1rbRFokedkmstON5GTkk5OurZqSUtJ5UFMMajmr4jI4dOgQ7du3r2ozqh2O7osQYo+Usruj8qpHoLDDL9CT8AEhdmk3TulAUv903L3cbJPcMVsT2frNUdzc9Zjyzfzz6yna9GhIgclMVmo+WWn5ZKflk52eT06GUXvgp+eTY+1NOGLXzyeJGNiU7iPDcPdUX02F4mqh/tsUl8Tgoadpe/vVEJ0HhtCxf2P0eh1rPtxP3P5klr24w6n2PLzd8A7w0JyKvzsFRjMn9yVjKZD889tpTu5PZuDEdrgZdJhNFho290ent5/OKuxZ5GQY7Z1NhhFpkXS6oQl1g30q7R4oFNcyyhEorhi99eHce3QLzh5PQ1okPnU98avrgU9dD3zqeOBT+MAP0B76Pv4e6A2l1yicO5lB2vkcotafIjUxm+9fj7LlefoYaB4Z5HTPAuDAljN0vSmUwMY+2qa9dCNmsyRiUAi+dT0r/2YoFDWYWukIzgQe5lxOEmpxZOVQr4kv/37j+gpNFjds7k/D5v607FKfbd+d4OS+C2Rd1FYp5WWbOPRX6b2GJXsW3gGao0k9m8Xh7UnsXhtXqs7e307TunsD2vdtTEj7umqCW6GgljqCVL+zvB71Ks9UtSHXEJX1QHVz19N/fBv6j29DbqaRvRvi0emFXc/CJ8ADbz93hz2LQtr3bcTe3+MROoG3nzvnT2dyPi4DgGO7z3Ns93na9GpIjxHN8fQ1qD0TilpNrXQEAMl556raBMUl8PJzp8+tLa+obuPWdUtpSF1MymbTsiMkHtOWxR7deY6jO7XvQVjnIEz5Zuo18SEsPIic9Hyy042Y8s00ahVAYCMf63yEEWNuAcEtAvAP8sSUZ7abp8jPKSAsvJ4aflLUKGqtIxDSAJgvWU5x7VA32IfbnuxKXpaJT2dsQecmsBRo8wxx+5MBOHPkIvv/SCivGRt6Nx3mgtJCgX9af7ftFUyXm5rh4W0gJ0Mb5qrf1M+m+aS4dnGlDPVzzz3H4sWLuXjxIllZWZXSZq11BD76IDQBVEVtw9PXwEMfDQQgOT6L+EOpePu7k5WWz6noFG0oqo4HPgHuCCGIP5xKTrrRNiwldIKzx9PIzynAzaCzpZuMZpLji/4xj+xM4sjOJLtzt+sdzKBJ7ZUzUFwxt9xyC//5z39K6RBVhFrrCAI9A1GOoPZSOKdRv5kf9Zv52dK7Dw9zqr7FrMV5KCkUmJdt4uLZbI5HnceYZ+bk3gvo3HTkZhgBOLwjicM7krjujta07tEQLz+DmrCu4VxNGWqA3r17V/o11FpH4KFTa8wVV45Or8NdX3qy2tPHQKNWdWjUqo6WMEnb3Wk2W1j38QHbENTWlcfYuvKYbXjqrhd7EdhIfScrxJwAF7VbdgTDqy1D7SpqrSMQVr29bFMOe8/vJbJBpMNyOaYcDDoDBr1aVaK4cvR6HcOnhXP2eBpnj6ex86eTALY5iq9e2glAy6716T6iOUEhvlVmq8J5rrYMtauotY6gkMSsRGb8MpHoydGl8pJzk7n1x1tp5NOIlbesrALrFNcSOp2gSZu6NGlTlw7XNcGYW8Cx3efYtfqkrcyJqAuciLrAkHs60KZncBVaWwMp583dVVxtGWpXoWSorWw4tcHuWErJi9teJD0/ncOph6vIKsW1ire/O3UaetNjZHMe+mggw6eFc9O/O9IgzB+A3z47SIFJrWqr7lxtGWpXoRyBlembppNpzLQdbzi9gc0Jm6vQIkVtQQhBi8j6tO7RkKH3dbSlm10cw1pRcYrLUEdERPDEE0+UKlMoQz1kyBDatWtnS585cybh4eF06tSJ/v37ExERwddff02nTp2IjIzk8OHDTJo0qVR7Tz31FCEhIeTk5BASEsKcOXMqfB21Uob6cFgWg3N78cDbf3E6CGbcp42QTewwkad6PIXJbGL0j6OJz4wHwKAzEDUxqrzmFYpKY9H0PzHmmfH0MXD7rO4E1FcxpctCyVA75nJlqGttjyDPUnojRkxyDADfHP2G+Mx4mvg2udpmKRQENtZWD+Vlm1j6/HaWz9nB9h9OVLFVimuZWjtZnJJfeg9B1PkoPt+3kqVHPgNgerfpzPxz5tU2TVHLufXJrpyLTee3zw6SdTGfi0k5XFx3igunM5EWiU6vw8vPQGinerTu3rCqzVVcA9RaR5Bt1tbt+rn7Anm29Df3vgxAqzqtGNh0YFWYpqjl6PU6Greuy+R5/Ug9m21bWhp/0H4i8siOJDZ/dZTAxj7c+kQXhBAUmMxkpuShd9PhH6SGlBTOUWsdQSF1E7N4NfRpnj31hl36hLYTEagdn4qqJbCRDw9+OJDMlDwunM7EzV2HlLD2Q201SV62icRjaXz44Ea8A9zJSTfa6g69rxOtujWoKtMVNYha7wgA6j32HpSY7J/5pZ6DfQ9VjUEKRTF0OkFAfS+7SeOJ/+1DSmK2zSEA5KQbETphC9izftEB8nPakpmaR0ZyHhnJuZgLLAy5p6PaxaywQzkCoG5+Dub8xug9zgNgzgsG6cbn2+Lwa3eJygpFFeAf5IV/kBcPLxiE2WQh+UwWXr4GfOt6YMw1s+zFHeRlm9i07Eipul+9tJMh93SgdfeGSvxOAdTiVUMlyYl9AuPFXgDknR1jlycl5BgLMJqNWKRa262oXugNOhqG+eMf5IVOr8PT18C//tuHbsNDadsrmO4jwxg8uT2hnerZ6vz22UHWLoimpi0frynMmTOH119/vdLbzcnJYeTIkbRr146OHTsya9asSmlX9QiAc15aAJP8pNG4pY3Akudhl18gTTz+zXZiPV+mgVcDlo1c5qgZhaLa4OHlRu/R9kF9WnVrwG+fHyT2nwuAFoMh5Uy20jWqYcyYMYOBAwdiNBoZPHgwv/zyC8OHD69Qm7W2R+BdtFCIM7718fNw48h/RxAzZwz759zEpD6hduW3mR8iKTuJ/cn7ueX7W8gx5QCaVtEn0Z9wJLV0F1yhqE64uesZ/kA4Dy8YRHALTakzP8dUxVbVfBYvXkznzp2JiIhg4sSJpfIXLVpEjx49iIiIYOzYseTkaM+OlStX0qlTJyIiIujfvz+gqZn27NmTyMhIOnfuzLFjx+za8vb2ZuBAbTWju7s7Xbt2JSGh4nL6tbZHcKaol4xZp+PXJ/rj4aYHwN/TwMujO7F4+ymHdeMy4ui1vBejW47mxxM/AvBO1Dt2wnWnM05zNvssvRr1ct1FKBRXSKFO2oYvDnHbjK74Bdb80JrhX4a7pF1HgpSFVKUMdVpaGqtXr+axxx6r4BXW4h5Blrfgrf6a9HRoPR8aBVz+mutCJ1DIkdQjZBozuXHljYz8fiT3/novv8b9CkCmMZNz2SpOsqJ64Gt98Gem5rH42W0c33O+ii2qmTgrQ3399dcTHh7OsmXLiInRFAwKZagXLVpke+D36dOHV199lddee41Tp07h5eX4uVRQUMBdd93Fo48+SosWLSp8HbWuRzA4O5dCLdF0Ny14iK+H3mHZz6f04FFtLw+ywJvcs7fj3XRxmW3fvvp2WtVpxbmcogf+4oOL+f3U7/wS9wsAXw77kq4Nu1b8QhSKCjBgQlvOHk8j66IWS3n9ogPodOG06FK/ii27csp7c3cVVSVDff/999O6dWumT59eKddR63oEIQUFRQdmbwB83B37w4HtGjD/urfIPz+UrGOzMWd1ICduGl2DrrOVyTl1r12d42nH7Y73XdhncwIAk9dNZvIvkzFZisZmk3OTMVuU5LDi6uHu5cakV/vSf3wbW9ovH0fz/RtRHPwrkbxsNXfgDFUhQz179mzS09N5++23K+06XOoIhBDDhBBHhBDHhRCl1jkJIQKEEKuFEPuEEDFCiKmutKckFrO2qcbL3XGPAGBY88GEud1CRNNAWjXwxZwbxtmj48k8NJ/MQ/Mx57Qi+4T9brTlI5Yzo/uMMtuMOh9FTHIMBZYC3trzFgO/GUjkEscR0hQKVyGEoMP1jblxagdbWuKxNDYuOcynT27h3MmMKrSuZnC1ZagTEhKYO3cuBw8epGvXrkRGRvLJJ59U+DpcJkMthNADR4EhaFHi/wbuklIeLFbmWSBASvm0EKI+cAQIllIaHbUJFZehbtskjmcitK5cuz8n8PK2xfj0v55mCxeWWVdKidkiafXcL2WW8Wuv+bkHOjzNf3r8C4u0MOK7EZzJOsMPo38g0DOQ/l/3L9fGDwZ/QP+Q8ssoFK7AmFvAj+/sxdPHjdMx2ltt+IAQux5DdUTJUDumOslQ9wSOSyljrQ/2FcDoEmUk4Ce0QTZfIBUowIWYig3n3dvPuS+5EAI3B4HK540JZ3gnLZxg5pE5ZJ94gre/1yaLdELHurHriJ4cTcs6LanrWZeNd24s9zwPb3jYyatQKCoXdy837pjVnVseiaTbMG3pdPSmBLZ8c5SMlNwqtk7halzpCJoA8cWOE6xpxXkfaA8kAtHAY1KW3rorhLhfCLFbCLH7woULFTIqR1d0yf1aBV1W3aeHFXXr/pw5gLt6NmPW8HY8f3MHsHhiMTagwCLJL3A83h/kFcS3t3xrO27g3YDvR31PA+8iYbDwL8NtP3vP72Xrma12kdMUClfTqpi09f4/Eljy3Hay0/Or0CKFq3HlqiFHU+klx6GGAnuBQUBL4DchxBYppd3gpJRyIbAQtKGhihhVEYGIaTe0oHeLQNo38sfToM0rhNbz4d/XNadzSAB3LNgOwIXMfELqejtso21gW/ZO3IteVzQvsXT4Um5adVOpshN/sd+cMrb1WOb0nVOBK1AoLk1QiC+jp0fyy4JojHnaS80XT//FtA8GoHfQM1bUfFz5V00AmhY7DkF78y/OVOA7qXEcOAm4VOatndFEv5xcmid1u+y6Qgi6NKtrcwLF6RFWtH740Nny3+CLOwGARr6NGB526S3iq46twiItbE/czq6zu5y0WqG4fELaBXLf2zcQPiDElrbqtT02ZVPFtYUrewR/A62FEM2BM8B4YEKJMqeBwcAWIURDoC0Q60KbcJOSBecucJexw6ULXyHf/5PAoHYNeHfDMVb8fZpl9/amcR1PUrKMNA3UegpSSqSE5Kx80nJNRHr+h9cmvYYQAiklf8T/wdrYtUSdjyI5N9nWdsTiCLtzhfmH8e2ob/HQ2+sjKRSVwfXjWhPU1JeNSw5z4XQmHz/6J8Me6ERwiwA8fQxVbZ6iknCZI5BSFggh/gOsB/TAZ1LKGCHENGv+AuAV4AshRDTaUNLTUsrkMhutRMb2bOWyttdGJ7E2eq3t+MY3/7R9ruttoIGfJ0fOle41/H7oHJ9N6YEQgsHNBjO42WBbXlnb5+My4ui+tDsPRjzIQ5EPVeJVKBTWJab9GrPjx1hyM4yYCyys+WA/zSOCGPFg56o2T1FJuHTAT0q5VkrZRkrZUko515q2wOoEkFImSilvklKGSyk7SSmXutKe4tzevVmlt3lPv+aXLHMxx+TQCQD8cfg8a6PPWnsL9l3wn2/7GdB6AENCh5Sq+9G+j+wmmtUGNUVlMnluX64f1xo3d+2RcXJfMlHrT3E6JqWKLaueuEqGGmDYsGFERETQsWNHpk2bVq4ekbOomZ9KZPqQ1nbHL4/u6FQ9t2LBQR5aFkXzZ9bS/Jm1pGYXbacI9Q9l36R9rL5tNW8OeJPoydHsm7SPAU0HOGxz2SElla2oPPQGHZ0HNuWOWT1sadu/P8Hq9/aRnZ5PRnIuFrOK1XE1+Oabb9i3bx8HDhzgwoULrFy5ssJtKkdQifh7Gtg2axA3tm/InzMHMKlPGMfnDmfOLR3Y9dxgtjw1kJs7N2L7M4M4OW8EcfNHEjd/JP+8UPoNH6DrK7/ZOQOdsP9z6YSOtwa85bDu/3b/j/Vx61XgEUWlEtjYh4ET29HphqKV4F88/RdLZm/n109iqtCyquNqylAD+Pv7A5rwnNFovKTWkTPUOtE5G+digLBKb7ZxHS8+mVy0ec9Nr2NKsSGj9yeUFpzz8zQQN38kg9/YxIkL2XZ5XV/5jeg5N5FwMZdPtpxkZOdgNhw6zz+n0xjfsykTe4faiW3N2zmP5YeXAzDjT3uZi/Ftx+Pp5omPwYdbW91KA+8GpZxLSZwR1VLULjr0awz9oHlEEKvf3WdLP/HPBU7uu0BwywC8fN2vul2H2rlmh3H7w2XHLq8qGeqhQ4eya9cuhg8fzu23317ha6y9jmD1o9Dnqk1JOMWGJwcAkJSeR+95G2zp4XN+tX1eFVUUhOKFH2N44ccYDr8yzLakdUb3GTZHUJIVR1bYPn+w9wNA2+T2w+gfSMlL4e09b7MxvvTu5x9H/6SGZdoAACAASURBVEiLOhWXulVcWzTrUI+HFwwiNTGbr17WZHrXfqS9lEx7fwB6t2t/wMFZGerZs2eTlpZGVlYWQ4cOBYpkqO+8807GjNHC4/bp04e5c+eSkJDAmDFjaN26dan2ANavX09eXh533303f/zxB0OGOB5VcJba6wgA8tKq2gKHBAd48tSwtvzfOueinrV7fh0HXx6Kt7sbBr2BfZP2YZZmxv40lpPpJ8utm5ybzHUrriu3zEvbX+LL4V86bb+idhHY2IeB/2rHzp9iycnQhjIX/GcTgY19uHFKB+o387sqdpT35u4qqkqGGsDT05NRo0bx448/VtgRXPsuuzwMlx+M5mrx0IBW/PxI0QN69sj2/P3cjayc1ocTr45g1YN97cpPX7HX9lkndBh0Bn669Sf2T9pP9ORodk7Yafu8bMQy7u98v9O2RJ2PIjbdpds7FDWcDtc1Zur/XUdQ06L4x6mJ2Xzz6t+cP3XtqphebRnqrKwszp49C2hzBGvXrrVTNL1SanePwKEKRvWhU5MA4uaPtEur76dtHOsWWpcdzwy2DSH9evAcp1KyCa3nY1e+8G3F21AkedG5fmc61+/Mw5EPc/+v99PYtzGPdn2UIC977aVfTv7CU5ufAmD0D6MZ13Ycz/Z69pLzCoray53P9GD/xgRO7k/mzJGLAKx+bx9T5vVDb7j2vjfFZaj1ej1dunThiy++sCtTKEMdGhpKeHg4mZna8vGZM2dy7NgxpJQMHjyYiIgI5s+fz9KlSzEYDAQHB/PCCy/YtZWdnc2oUaPIz8/HbDYzaNAgpk2bVuHrcJkMtauoqAz1oIbH6RKoedSsnp8R/8TsS8pQV2e+2R3PU9/avzU0DvBk2zODy6hxeTjayFYVkaAUNY+tK4+xb4OmO3ndHa2JGNz0EjUuHyVD7ZjqJEOtuAqMimhcKi0xPY+wWWuYtmRPhduPnhzNzS1utkvLNmWXUVqhKKL4g3/rymNsXHqYpNj0KrRIURbKEdRwPA164uaPpG3D0hNy62KS2HbcsWKHyWxhyfY43t1wjIeXRxE2a43t59eYJJKzimSH510/j6i7i+Ygei/vTa9lvfjswGeVfj2Kawe/QE+G3tfJdnxwayKr/m8PH0z7g50/qTmn6kQtnyO4dlj/uLYhJfZCFoPeKNI2mvDJTmJeGoqPhxtSSvILLHy69ST/W1/2iqT7y+hJ+BXraeYU5PDWnrcYFjaMxr6leyUKBUDLrvUZN7sn0ZsSOLi1SHx499o4etzcHJ2ues/T1RaUI7jGaFHfl7j5Iwl/cT2Z+Vqwt44vrr9kPU+DjjxT+RIBOafvwbuZfS9g6KqheOo9eXvg2/Rt3FdtPlPYIYQgKMSXgf9qR/9xbUhJzGLlvMuf41O4ltrtCIzX7lj37udvpO3sdWXmP3BDCx6/sQ0FFomPu972ALdYJMfOZ/Hc99HsPnXRro45uw2Zh+YD4NlkKQb/AwDkmfOY9ru2cqGD53jCfW7jvv7NOJUVS5fg9hh0RXLFBZYC3HS1+2tXW9EbdDQI9S9KkJLqvnKvtlC7/yP/eruqLXAZHm56js0dzuvrj/DxZm089oH+LfD1cOO+/i0cBtcB0OkEbYP9+LbEPoU9p1IZ+9F223Fe4jik2Qf3ujvtyh3MW8HBvBV8/X359r054E2HKqqK2sPe3+PpOjS0qs1QUNsdQXoCUK+qrXAZBr2OZ0a055kRFV9e1y00kLj5IykwW9iXkM7Yj7aRn3Qb+Um3ARL3epvwaHDpIahCntj0BANDhvDu4DcrbJuiZiF0AmmRbP/+BOEDQzC4O34pqcnMmTMHX19fZsyYcenCV8ioUaOIjY3lwIEDFW6rdjsCxWXjptfRLbRuqY1u5zJu5GzOVCauH8+wZrdx9GgE++ONuNfbhNBnY8lvhDm3Kd6hi2x1Nib8poaKaiF3PtuDr/+rhVqN259M6+4Nq9iimsd3332Hr6/vpQs6ifoPVFQKDf09aejfsWiz2UA4k5aLxTKaHKOZDzYeZ2y3ELo2e5Bu78/Eo762I7rLki4ADK4/jTmD7qaOZ52qugTFVSIoxJfmEUGc3JfMr5/EUDfYh6CQynuoXW0WL17M66+/jhCCzp07s2TJErv8RYsWsXDhQoxGI61atWLJkiV4e3uzcuVKXnrpJfR6PQEBAWzevJmYmBimTp2K0WjEYrGwatWqUsJzWVlZvPnmmyxcuJA777yzUq6hdjsCnYq56kqa1CnScnr3ri62z7/f+zwjf9xgV3bDhQVs+HpBmW21qtOG2b2fpVWdVgR4BFS+sYqryvXj2nByn7bHZfW7e5n6f+ULHzrDB9P+qHAbjnh4gWPRN6gaGernn3+eJ598Em9v71J5V0rt3lDmUXPfQmoyzerUY3aHNWSffMTpOsfTjjJl3RSuW3Ed4V+Gk5Wf50ILFa7GL9CTbsO1ieKcDCNfPL2VAmPNC6/qrAz19ddfT3h4OMuWLSMmRgvgUyhDvWjRItsDv0+fPrz66qu89tprnDp1Ci8ve2HMvXv3cvz4cW677bZKvY5a2SPYdqEZKfneDKhjAnKr2pxaybgezRjX437gfkxmE9/9E8cLW+fhXvdvTBmdEfpM3HzKltDus0ILmZhz6t883Hs4vVvWo2/LoDLLK6ofPW9pwdnj6SQeSyM73cjHj/7JA+/egNsVTh6X9+buKq62DPX27dvZs2cPYWFhFBQUcP78eQYMGGBr80qplY5ge7L2JhIZ4gVUfMZdUTEMegPjurdmXHf7zWr5BWY2H00mO9/Ei+s20KlxHfYz266Md+invLe1Du/+Ub9Uu6se7Eu30LoutV1x5eh0gluf6MJvnx3k2N/nANi1+iR9x7aqYsucZ/Dgwdx22208/vjj1KtXj9TU1FK9gpIy1E2aaGE+C2Woe/XqxerVq4mPjyc9Pd0mQx0bG8v+/fvtHMGDDz7Igw8+CEBcXBw333xzhZ0A1EZHMHQuLP4cAJ3Bs4qNUZSHh5ueIR20FSW3dpkMgJSjeGHry/wQ+62tnG/LNwDIPPo8mItkuMd+tM32+Y8nb6BFfTUUWN0QQjDg7rY2R3A86jx9xrSsMTvUr7YMtauodTLUA6c8wMYvPgbg7shALn75d42Woa6tpOSmMPansaTkpdil+9OOxMNTcPS1ntCrGXqhbZi7qWNDdEJQ11uLratXmjdVyqFtifyx+DAA7XoHM3hKB+fqKRlqh1yuDHWt6xGYTUbb5+L/+lJKLBkZnH/9DXIPHKD5N18jDGpVUXWlnlc9No3bRFJ2EkO+LdqhnMFhfNvNIsA9iHY+AzmRpONkrLZiafnO07Zys39wPCTYPMiH/wxsRWSzOrRUPYirRljnovmdwzuSuBCfxR3Pdkevr93rWa4Wtc4RmMpYbZL08sukfVUU3N2YkIBH8+ZXyyzFFRLsE0z05Gg2nNrA9E3TbenpxmR2GleCB4SE/8mZwxNBuiHN1lUY0t1heyeTs3ly5T67tGdHtOOefs1xUw8ll+Hl687DCwbx/RtRJB5LI+VMFif3JtOqW4OqNq1WUAsdQX6ptOzNW0ql6dwdPygU1ZPBoYP5666/eOGvF9hwegP1vepzIfcCAOkFSfi2+l+Zdb8bto0vth9hxYFNmHPCwGK/ZO/VtYd5de1hu7R2wX48Org1wzsF15jx7JrAqMciWfCfTYAWzCa4hT++ddVcnqtRjqAMzsyYSejyZeqfvAbh7+7P2wOLhATjM+IZ8f2IS9Ybs04T2PO2BtR6oPMDNNON4dMtsexLcBxR63BSJg8ti7Idzx7ZnsHtG9I8yMdheYVz6N10dBseyp5fTpGdls+Xz2zjvrf74+5Z6x5VV5Vad3cLnHQEuf/8Q96+fXhFRrrYIoWraOrflN3/2o1e6Dly8QghviGk5qXS0Lshsemx3LXmLof1Pt7/MfAxYU3C2H//Mnzd/Ii/mMOBMxks2hLL3vi0UnX+u+YQ/11zyC5tTNcmzL01HK9rUFTNlfQY0Rwk7Fl3CoBF0zfz0EcD1UuZC6l1jsBkdM4RAKSv/lk5ghqOh94DgI71OgLY5Ck6BXUienI0u87uIi4jjkHNBpGUnWTnHOIy4uj3VT98Db78Oe5PQus1YmTnRrb8A2fSufm9rWWe+7uoM3wXdQaAf1/XnJ7NAxnaMbjSr/FaQ2/Q0fvWlpw5epGk2AwA/lhymMGT1OogV1HrHEHBZUgTZG0pPXeguLbo2agnPRv1BCDIK4joydHEJMcwfs14W5ksUxbdlnYj1D+UVaNW2ZxLpyYBNhXWC5n5LNt5irjkbM6k5fJ3nH1Qn0+3nuTTrdpO6ZPzRqi3Wye47cmufP7UX+Rlm0g/n1PV5lwWrpShHjBgAGfPnrXJT/z66680aFCxSXWXOgIhxDDgHUAPfCKlnO+gzADgbcAAJEspb3ClTcXnCC61h0IajeXmK65NOgZ1ZP+k/aw/tZ6Zf860pZ/KOEX3pfbLsHdO2Im3wZv6fh5Mv7GNXZ7JbKH1c7+Uar/5M2vZ98JNBHir5cnlodPrGD4tnO/fiOLs8XQyU/PwC1QTxwDLli2je3eHWwKuCJethxNC6IEPgOFAB+AuIUSHEmXqAB8Co6SUHYE7XGVPIcWXj5pzTa4+naKGIoRgWNgwoidH89vtv5VZrtfyXuw6u8thnkGvI27+SOLmj+TQy8Ps8iJe/pWwWWs4mJhRqXZfa9RpWKSwef5U9bxXixcvpnPnzkRERDBx4sRS+YsWLaJHjx5EREQwduxYcnK03s3KlSvp1KkTERER9O/fH9DUTHv27ElkZCSdO3fm2LFjV+UaXNkj6Akcl1LGAgghVgCjgYPFykwAvpNSngaQUp53oT2A/WTx1tPZtHHT41lQpHro1iCI0CVLOTF0mKPqilpI4V4FgLj0OGLTY3ls42O2/H//+m/chBut67bmUOoh/t3p39zZ9k4a+TSyDQF5uevZP+cmOs/51a7tEe9qw4+jIxvz1p2R6NQOZzu8/d0JC69HXHTKJcu+Me5ml9jw5Nc/l5lXFTLUAFOnTkWv1zN27Fhmz55d4aFGV+6QaQLEFztOsKYVpw1QVwixSQixRwgxyVFDQoj7hRC7hRC7L1y4UCGjivcIEvIkMSH2ipVBLc8gLmMPQX5sLIlPP43pzJkK2aWoGYQFhDGo2SCiJ0fzQp8iHZgCWcChVG3V0KcHPmXoqqF0XtyZ8C/DCf8ynHxzPv6eBlsP4d7r7Dcr/rg3kRbPruWb3fEo7NFZN/Kt+/gA6z+pXiKRV1uGGrRhoejoaLZs2cKWLVtKBcK5EpzqEQghPICxQFjxOlLKl8ur5iCt5KC8G9ANGAx4AduFEDuklEftKkm5EFgImtaQMzaXRcl9BPl160CcJngVNuQCnoEmCi6jvdgR2mShzseH4KskEKWoHtzR5g5aBLRgyroplyxbcm4B4JWJz/F/a05jdD+KMeUGpCmQp77dz1Pf7gfU8tNCWnVvwJmjF8nPKeD47vP0uiXHbsiokPLe3F3F1ZahBmzqpX5+fkyYMIFdu3YxaZLDd2incXZo6EcgHdgDOLv+MgFoWuw4BEh0UCZZSpkNZAshNgMRwFFcRElH4OFWNGHnVc9+zqAgKQnTufPk/hPF+TfeJOTdd/AsJuSUd6TITF0lRgtS1By6NexWFJ7TSuEihIMpB+1WH5Xk/3bPhYbgDrjX3QmAMbUPBRkRmHPD7JafAjQL9OZ0qja+3LN5IBN7hzKgbX38PK/tSefW3RvSqmsDPnxoIwC/fBzN+Nk9EdVgGO1qy1AXFBSQlpZGUFAQJpOJn3/+mRtvvLHC1+GsIwiRUl7uoPnfQGshRHPgDDAebU6gOD8C7wsh3ND+H3oBb13meS6LkhvKPLy0B7hwszgsf3ryZIxxcQCcvG0MbXbvRu/rQ8a6daQs+sRWTufn7xqDFTWOwjfEjkFaDOfErETm7ZyHl8GLEN8QjqcdZ2P8Rod13QO34x64HQCLsQ65Z+7Gkqe9TxU6AYBdJ1PZdbJoPLpX80Du6tmMW7uUHH29NhA6wajpkfz09l5SE7PZ/UscPUZWvRbY1Zahzs/PZ+jQoZhMJsxmMzfeeKNt/qEiOCVDLYRYCLwnpYy+ZGH7eiPQlobqgc+klHOFENMApJQLrGVmAlMBC9oS07fLag8qLkNdkk7XD6Tr1gUENEnB3U8bpzM9cJjjAx1HO2r031fQeXtz5okn7dLrP/44QQ/cf9l2KRSnM04z8vuRTpeXZi+yjs0C6eEw/5XRHbmtawi+HtfeNqEFj2zCbNJe2gbc3RZdUIaSoXaAq2SorwOmCCFOog0NCUBKKTuXV0lKuRZYWyJtQYnj/wFlK4K5mANbNuJVz5v+ftYFS+7lSw/rAwO58N77V8EyRW2hmX8z2/CSRVrYGL+R6Runl1le6HPxa/ciAC19u3B4753ky3R0Htp3+PnV+Tz/YwyT+oTy5JC2CB0cSEhHrxN0Dwus0bEXxs7sxjev/g3ApmVHGPhoo0vUUDiDs45guEutqGL+TmlK/wZxtmNZxpItAEt2NvmHijRlPDq0J//goTLLKxSXg07oGNxssM0x5JhyOJlxkhxTDpviN7H44GK78iey/sHQ6h8czRJ8nwZfvTcKU3qXUoqqTw5pwyODW7vqMlxG/WZ+DLmnA799pq1Cz882YTFbbCuLFFeGU3dPSnkKqAPcYv2pY027JnGrV6/MvMSZT9kd+/bt62pzFLUYb4M3Het1pEdwD2b2mEn05Gj+vvtvpnaa6lR9z+Cf8Gv7En7tZ4GuaH7sjd+OEjZrDb1f3cB3UQmcz8izTXIbCyzkmcp+GapqWvdoaPtszDOTkph1SZWA2sSV3Atnl48+BtwHfGdNWiqEWCilfO+yz1hNWXW6I7c1jUFnzEJ3YR8hH35IwkMPAWAICcG9RXO7uAV1xo8j6IEHuLj8q6oyWVFL8XTz5IluT/BEtyfINmVjlmb83bXFCmaLmeWHlxN1LorfT/9uV8+v7YtE+NzFtpg6WPIbgsVAUkYeT3yzz9FpaB7kw2+P9692AXmEENzz+nXEbEkkM/UCGfXS8c/zwsPr2l495QxSSlJSUvD0vDwpDmeHhv4N9LIu80QI8RqwHbhmHEFcdiCp+d4EeebAZ0PRDyh6wPtcfx0Za+01Yxo+/TQ6B5s9FIqriY/BPv6BXqdnYoeJTOygSR2YLWb6rehHtikbgH3ZX+ET5rgtaTFgzg3FktcIU3pXTiY3pJVVK+mB/i14aECraqOP5OXrTvfhYXw8/SSQREb2RTx9q4dtVY2npychISGXVcdZRyCA4n1FM443jNVoDLpil7hyMqDtFvTp08cujKVH+/bKCShqBHqdnh0TdjDjzxmsj1tfblmhM+Hmcxx8juNer6j3K80eLL+Qz/KVkJ98Ax7ZN1Lfuw6jIhvzQP+WVbrhrfMNzYj65RSBjS2Mmx2uJDquEGcdwefATiHE99bjW4FPXWNS1VF8g6CloOjAp2dP6k+fzoW3tZWtgXeX3A6hkRsdTcG5c/hVwgYPhaIyef2G13n9htdtxxZpYeWRlRxIOcDRi0cJ9g7mj/g/HNYV+qK5BY+gPyHoTy4An56BT78C48We5F8YSudGjfl8Sg/c3XRXbZNb3WBtH1BqYjYfPbSRe/53HV5+Kszs5eKUI5BSvimE2IS2jFQAU6WU/7jSsCrBsy6YkwDQGYo2mOm3zCFo2tt4d+9G3tGjBIwdW6qqJTODuDvuBKD1tr9wc6A5olBUF3RCx7h24xjHOIf5eQV5zN81H5PFRBPfJny076My23Kvuwv3urs4CfR+/zZMab0AWH5fL/q2DCqzXmXQokt9zhy5yOEd2v/tynm7GTS5PSFt67r0vNca5W4oE0L4SykzhBAOn2pSytJSey6msjeUFee+9z/D/7s74OxepIQMw0g8U9fh4W+GOY5j155/401SFi1CHxiI2ao8GLp8Gd5du162jQpFdafAUsDvp3/nRNoJlh5cSpYpy2G5zEP2oUc2zhjg0njO3/1vD2dPFP2P3vJIBE3a1UVfzSa6q5KKbChbDtyMpjFU3GMI63GLSrGwOvHAn/BGe0RmIgF1jkGBc8vozMXkZ09NuJt2hw6qKFSKaw43nRvDwjS1mYcjH7alb07YzMMbio792s8CwJzXiNzT9zLw9U0A+Hu6kZFXQAM/D/Q6wYr7exNar+IOou/YVqx+dy/GPO3/dfV72kqo7iPC6HJTM9w9r71d1pVJuXdHSnmz9XfVi3pcTXzqQWYiJLtM+06huKboH9KfqIlRdF1i3xPWe57Ft80rAJjSupKRdCvgzvlMbd7hhv9t0srpBHtfGHLFcwvBLQK47+0bOHUghZ/fL1oOu3ttHLvXxnH7rO40DFN6YGXh7D6CfsBeKWW2EOJfQFfg7cKAMtc8nnWq2gKFotpj0BnYP2k/P8f+zLJDy4hJibHPrxOFoU6U7TgvaRSmtJ4g3TBbJOHWoD09mwdiMluY0jeM0ZGXJ6IX2qkeDy8YRPqFXDavOMrpGC2gzbmT6coRlIOzA2gfATlCiAjgKeAUUPFoCDUF90t3Xd1btKDFL2svWU6huJYRQnBLy1tYcfMKoidHs2/SPsKDwh2W9Qz+Cb92s/FrPwvhVjS+v+tkKv+cTuOxFXsJm7WGsFlrWLP/LBl5Jqd3zQbU9+KWRyIIH6Ctp89KdVY9v3bi7MBZgZRSCiFGA+9IKT8VQkx2pWFVSoHzXxqvLl3Q+frSYMYMPJrXrhE0heJS6ISO5SOXA9rmtujkaH44/gOrjq2yK+fbeh4AgZkPQm57TqXk2OU/vLyoJ+Fp0LH/xaG4u136PbZeE+0l7p/fTtNzVHPcDLU7yE9ZONsjyBRCPAP8C1hjDUx/7W7jKzk3kHEGFg6EvNLBs/0GDaTN37vwGzTwKhmnUNRM9Do9kQ0imdN3DtGTo9k6fitdG9jPKaT6fURqg0fZNrsL66f3d9hOnslCm9m/EDZrDUnpeQ7LFNK+b5E6afQmFU62LJx1BOPQ5Kf/LaVMQos9XGXS0VcNv2ISt4lREL3SYTG1OkihuHwCPAL4cviX7J+0nyGhQ+zyhq4ayu3r+/Lh/cIW53nP7NIbNXvP20DYrDW0f34d6w4klcrX6XV4+WnvrNtWHScnw+iai6nhOKs+miSlfFNKucV6fFpKufhS9aorTRv58Xi7LTzebgvTl36PX1D9Mgr2sj9e84TrjVMoahlCCN4c8CZLhpeednx6y9OEfxlO+Jfh5MjzxM0fyaYZA0qVyzWZmbZ0j21OIWzWGh5cugeLRXLDXW1t5bZ+o1YCOqJcRyCE2Gr9nSmEyCj2kymEKD1OUkPQUYBOgE6A3uBghKvLRDD4QJ//lM4z5pROc0DGTz9V0EqFonYR2SCS6MnRRE+O5s42d5bKH/HdCMK/DOeWNb355EEPNjzdgSeGBdOrhePVQL8cSKLFs2sZ/PXfnNdpSgHHdp9HWpRkdUkutY/gOutvv6tjztVBZykov8Co92Dkm5B9vnRe1GJoOQh8gsC7bBmJxKdnETB6dAUtVShqJ8/3eZ7n+zzP2ti1vL/3feIz4+3yH9/0eNGBB/i1h/B6EXTyHkd+VjO+3JZQlC9go5eJcdlaaM8PH9rIMt98Cuoa+GJqTzo0VstKnd1H0BuIkVJmWo99gY5Syp2uNM5VCGm6RAEBbu7gaKnauqe136H9YKpaLqpQuJIRLUYwosUIjGYj/5z/h3t/vbfMstEp+4hO0TaT+bWHr4etxVLgz6trD7HteIpd2buzPPhW5jPi3SKV1VERjfnvbZ3wv0qCedUJZ5ePfoS2iayQHAdpNQad5RKOoJCAEOg6GRq0h3Wz7PNO/aX9ltJOttS7d29yduzArUGDSrJWoVC4693p1aiXLYRnIRZpYXHMYt7Y80apOuPWjQCgUd1GxM5fR0paPj8sPYgxRtuzMCDXwElD0VLxn/Yl8tO+RAAeHdya3s0D6R4W6NQy1ZqO0/EIZLGdHFJKixCixop3CHdNuvaSO4aFgFHvap/TTsOOD+3z5wRov+/6Gtpq+iuNX53L8UGDQa/WKysUrkYndEzpNIUpnaYAmmOIWBxhV+Zs9lkiFkcwsOlA5j4wl6N/pLDjh1iCLDq2Pz2IL3ecYsGfJ+zqvLvhGO8WO35pVEf+1TsU/TUa78BZVxcrhHhUCGGw/jwGxLrSMFeiaxIJvabBvRucr1TOfAA7PnC6mfSffiJr82bnz6tQKJxGJ3S2Ceffbv/NLm9j/Eb6ftWXuCZFCvoN/DyYNbwdcfNHsn/OTfRp4The+Ys/xdDy2bWEzVrDYyv+qdYxna8EZx3BNKAvcAZIAHoB97vKKFcj3L1h+GsQ1Mr5SlkOJo4LySy9ftkRGb/9RuJTT5P43HPOn1ehUFwRwT7BRE+O5qkeT+GuKwpWMyfqBcxCWzCy9/xeW7q/p4Gv7u9t27dw4tURtKxfWl7mx72JtHt+HWGz1pCZ5+QwczXH2cA054HxLrblqqG7kmGbvo9CUrQ2nHTUPn6xQ5VS60haQWoqibNmYcnOIf/wYS0rt/zdkAqFovIojOEcdS6KyevslXHuWX8PFp2ZedfPY3jYcPS6omeDXifY8OQAAPJMZuatPcSX20/Z1S8UyvPzcKNfqyDG92zKgLY1b37QqR6BEKKNEGKDEOKA9bizEGK2a01zHUJ3BZM/dZrCPetg4LNFaf0eK7N4QVISsYN6cqxvP7I3byF3zx4s2VoAcUtWFofatSdt1aoy6wPk/P03Z2Y+hel8Ob0RhULhFF0bdmXbXdsI8w9DJ7VnQFC2Jkr3zJZniFwSyfifx2ORllJ1PQ16Xhrdibj5I9n3wk2l8jPzC1gXk8SUz/8mbNYaDifVrG1Wzj4RFwHPACYAKeV+anAP4Yp6BIU07Aidx8PNb0G/6UXp5w9BXrrdktP8otquIAAAIABJREFUxMxymzr73Gyk0X7Lu8VoJDc6mqwtWzn9wDQyVq8ma+OmK7dXoVDY8HP3Y/VtqxHWR9+YA0/gXuBly49JiSFicQRfHPgCs8XxPECAt4G4+SM5/MowVk7rw109m9Iu2H6r1bC3t/D0t/tddyGVjLMrf7yllLtKaOpcYldW9UV3JT0CW2U9jPlY+2zKLUr/sLf2O1sPNLQlZ3oYcLNY8DKZqXPHHaSttNcrSvnsM+rdey9YLMRNuJu8AwdKn7OML6RCobgyug4LJWqdNsyz4eZNFPjnMPCbIuHIN/a8YVuSGuQVRF5BHgObDuSmsJvo07gPHnoPPA16eoQF0iNMW0gipeSuRTvYEatFK/x6dzzrYpLY+8KQaq9H5uwTMVkI0RJruEohxO3AWZdZ5WKErpKWdurK96MZnu5sbdOU/TdeR8hHHxL88ks0fuN1vHv3tpW58PY7HO4UzuHOEaWdQDX/8igUNZU+t7bEL9ATgI1LDxPkFcTfd//NiOYjSpVNzk0my5TF6tjVPPLHI3Rf2p1Xd75Kvtlerl4IwYr7+3DgpaG2tPRcE82fWVvtVxk56wgeBj4G2gkhzgDT0VYS1Ugq1CMojt4APvaCdSeMdUkLcEfvZeZE+wCkTpCTlYnfwIEIIQgYOZKQt9+6ZNPuLVvif8vNAJz7v/9xaspULHl5TgfmUCgU5eMbqElOJB5LY8UrO3HXefBa/9dsy08X3bSIhyMfZkDTAaXqfnX4K7ov7W4TxLuYd7GoXQ83jv53uF35ds+vI+KlX8nOr54DKZccGhJC6IDuUsobhRA+gK5QaqKmIipzs9eUNfD5CMhJ5nhmIGuS2kMYtA+2cDZJczi5GemsX/AuN977IHo3A/o6dWjyzjucecx+srnhC88TOGECBamp6AMCSHpZi/X6/+2dd3hUVdrAf2daMukJKaSQhEDoHaSjSFFAxbI2sDds2N21fYru6qrrFhuKrgV1RVCxYENRQSnSews1oSSE9GRSp5zvjzukk0ySmRTm/J7nPtx7T7nvScJ97znnLbK0lJK1a0kZNBj/sWOJuPceTEnd0Ae0POm3QuGtXPrgEFZ8ksLulenkHC/myM4cEgeEV5aPjB7JyOiRNdqkFqRy0VcX1enr7EVa7gQfvQ8ze81kXNw4Dj0/lSe+3MUn67WMvgWlVvrO+bGyzUOTezBjRDzhAT6eGF6TaPTTWErpAGY7z4s7uhIAN84IACJ6wl8OIufks6K8ao1xz4maz9i5/CeO791TeR04cQLx8+fT9asvCb/rThI//5ywmTMBMISFIfR67IUF1KZ41SpSr7yKzOf/7r4xKBReiNAJxs+sClGdtjOngdoaicGJlTOGuRPrOpKW28t5f9f73PzjzQz8cCDWTh+w/ekJ9fb1r2X7GPbsz5VhLdoSV9+Iy4QQDwshugghwk4djTUSQkwRQqQIIQ4IIR5toN5ZQgi7c+/B47TIaqgaRbnZSIcD6XCQtmMrBVmNmXlWLesIgwH/kSPw7dWLiHvvxdyvb53awRddhCkxEUN0dJ2ygsVfULJpU0uHoFB4NUIILr5/EAA7fz/OpqWpLrc9O+5sdtywgzUz1nDP4HsI8akbsubH1B8Zs2g4gb0f5bZLtvLN7LEAdPKvcnC795MtJD76HeW2tttHcNVq6Ga0t9hdte4nna6BM53lXGAymjfyBiHEEinl7nrqvQj8WLcXz3A6P4KyYgtBEa45g6z/+nNWLpjPhJtuJ2N/CntWrahT56LY3Xx7vDeS5m36Bk6YQOAE7WvCUVrKvtFjkKVVlkr5i7/Ab+jQZvWtUCg0YnuGVp6v/eoQvv5G+o6Ldbl9oCmQWQNmMWtAVbCFrJIsJnxWcyawMGUhC1MW8uptzzG922SW7c7ktg83Vpb3/L+lfHzrCMZ0D6e1cXVG0Aftpb4N2Aq8BtT9hK3JcOCAlPKQlLICWAjUF6D/HmAx0GpeU7rTWA199Mi97F+/ptH2UkpWLpgPwK/vv1VDCXSKiyd5xGhmdV9Hj6Acru26hSBjyz2JdWYziR//j/j336u859a9DoXCSxFCcM0zVXsBW5YdwdHC5DURfhGVS0hfTP+iRtkTq56g/wf9eW7HVbx3lxmocmC75p11vLOy9cO4uaoIPgB6A6+iKYHeznsNEQtUzyZxzHmvEiFELHApMK+hjoQQs4QQG4UQG7OyslwUuYH+GtgjWLWwbrq82hxP2V3v/b7nTOTGf73B9AcfJzBQc1KJ9C0m2A2KAMC3Tx/8R40i6kmnU7dBKQKFwh2ERPlxx9zx+PgZKDhZysqF7ktpmRyazPbrt/Pg0JqpbnPKcrhv+X0E9n6c0D5PIQz5ADz73Z76uvEoriqCnlLKW6WUy53HLKBnI23qWw+prWZfBh6RUja4OCalfFtKOUxKOSwi4jT5hZtAQ4qgU2yXRtvv+HlpnXuJg4Yy8eY7q27M3gTxo2tWyt4PrwyERdfVn/TGVZz+BWXbd+CotlSkUCiaj16vwzdAS0qTm1Hs1r6FENzU7yZ23LCDD6d+yOSEyTXKbbKCgOQXCOz9KHq/A0z81wq3Pr8xXN0j2CKEGCmlXAsghBgBrG6kzTGg+ls1Dqi9PT4MWOj0ugsHpgkhbFLKr1yUq1nU3iyWjqqpWWjM6dcGpcPBqoUfsnvl8hr3J916FwMn13JECYiAm3+A52Kq7n33IPgXQF4qFGdrdZqDc9patmsXKYOHIPz8kGXarCPywQc0T2WFQtFkzpnRkyWvbCV9fz5z7/gVgFGXdWPIeQlue8bgyMEMjhwMwKGCQ1zy1SXIat/IfgnvcBLo9epgll37Jl3CPG8m7uqMYASwRgiRKoRIBf4AzhFC7BBCnC6gxgYgWQjRVQhhQotNVCOju5Syq5QyUUqZCHwO3OVpJQB1zUctuVVmYz5+p/+hr1+ymPVffw5AQFgnonv04qpnXqyrBE6DxWZCSsgq82PjY+fgmBMMb46FH58AW4UW6rrcAgXHoOD4aWcN+sCAGteypAQcDnA4yHrl1XrbKBSKxjnlbVydP744SHF+eT21W05ScBLbb9jOjht28OyYZ2uUGYO3MO2bkdy3+GePPLs6rs4IpjS1YymlTQgxG80aSA+8J6XcJYS4w1ne4L6AJ2lOiIkNSxaz6pOqbZEJN95O8ojRDbRwYq2aYv6Q3ouTZQHsLoik1G6is7mIuMwdkLkD/ni9btuJc2Dcg1pMo31LIWk8mEMJPO88Oh06TM5bb9VtY+iwieMUijYnJMqPmU+PIC+jhOKCcn537hXMf3Q1d71xLsKDGcou7n4xU7tO5YX1L/DZvqqYZL9aHuD2j97nreuGeezZruYjSGu8Vr3tvge+r3WvXgUgpbyxOc9oDs1xKPv94/drXCcNHe5y27yKquiGm3LjKs9tDk2O4yWBHLJ0YlR4GgZdtVnAL89oR3Xm5KMzm4l84H7C77id8n378OnZk4rUVA5fcin60EbSbyoUigYJ7exPaGdtZcA/2Icf3tLyJL9x13IGT45n9J+akNCqiZj0Jp4a9RRPjXqKf66bxwd7Nae13/L/zbPf/h//d2Efjzz3zM/KXA9NdSgrysmucX3DS6+jd/XL+/L30Pn41VuUL8I5WebPwrRBrM/pwrGS4Mb7S1sDhVq8P53ZjHngQHS+vuh12tTVlp5B6tUzOPHsc1hWriT/iy+RtvYZ30ShaO8kDgwnJKrq/++uVa3nBfzwiKpwbsagHczfvhibvW6uBHfglYqgqYlpqvsWXPrIHMLjE11v3O9PTLrnyRq3zEHaC/+XI9F8dHhI5X3b1Z/C0wWU3bKKI9YYpNT2hWuYNM+fBgtnaOdFJ2DuCHg6GP5bFd6idOtW8v73P47eNouMxx9nb7/+nHzlFYp++QVpb99REBWK9oROp/kYzHx6BAAVpTYKslrPUu+PGX9UnptjP2V5SsvN5+vDKxVB7RmBwdRw0Kf96zRFMO2eh0kaclaTn9d18DBm/O0l9AYDI/90NQ776b/QSwoLWPDvuXx2oBt7Jn7Je7mX85+943j3wDDyK5wbWelbtJf/v3pClpb+0uBnJyDm9P4KOW/O49jds9nbtx+pV13Nnl692dOrN8Vr12H57TdKtmypYT2lUCiq8A+pekcc3uaZl3F9BJgC+OSCTyqv/++rrQ3Ubj5eubNYe0Zw62vv8NNbr3Jo84bKe2UWC3tWryC0cwzH9uxEbzCQNMT1fYHaxPTozez5n2EwGsk9fox9a1cBMO3eP7Nn5XIOb9nI1y/9rUabpW++XPlyzrea+Sl7IOeHbyXYVNeCQUyaQxeh7SeU5RvI3ByC77RZ5M6fX6du6bZtledHbryx8lwXEID/mDEETZ1CwNlno/Orf0lLofA2TL4G+o6LYdfKdLYuO0LysKgaysGT9O1UFcShJPZhyqyT8DW615nUKxVB7c1i/5BQOsXF11AEi//+JCcO7q+8ThgwGJ8WvhgNRs1Z5ZzrbqbHyLEkjxiFTqfnj88+rrd+7S/0owUmVgeNYZpJs2/mnEdg/GNVCWy6T4QPpuMbkk/ChGy4bjJhN1xPRWoq9sKiqrDXOp1mbloLh8VC0Y8/UvSjFvYp4v776XTLzQin3AqFNzPqsu4c35dPfmYJ8x9dzaUPDSYmObTxhi1ECEHP0J6k5KUQaAp0uxIAL10acsV8tLoSAEgeMcZtzw8Kj6TnqLGVMY+ikpJrlMf07FN5/9rnXyYisSq2X55PV3gkDebkw7mP18xiFj0Qbv+t6vq/52I8uAj/ESMIOv88eu/dQ7efl9Fzy2Z6bNxI/Pvv0XPrFnz71h82Kuvll9nbfwAnX37ZTSNXKDouPmYD591S9X/ly39toaKsdQwxPrngE54Z/QxPjXrKI/1754xAf3r9t3LB/MqActXpNmyEx+QZN/NGEgYMpqK0hKQhw6koLWH/+jUMnXYJvgEBTLjpdr54/mmsZaX4h4SBuQET0dBEiOoHmc60lz/P0Y4nMsHoiynOab7q44P/qFEAJH72KQhByfoNWI8eIeP/am5u58x7i9CrrsJYTzhshcKbiIgPZOwVyaz6TPtQXPTcBq7960iP5yQ26o1clnyZx/pXMwIX6Dn6bMwBgR6SBoLCI+g3fhJDpk4nJKozkYlJjLnyWnwDNA/iuF59mTpbC1h1cONaVnz4TsMd3ra87r3norQN5i0fw47PYc+3sHMx2MoROh1CCPxHDCfk8svpvXcPcfPerNH8wLkT2DdqNEduuZWSzZvdMm6FoiPSd1xMZUyiwqxSMg8XtrFELccrFUFDM4L6uPC+v3hIkuax6buvsNusp69gMGkzgLvX1y37+i5YfAssugY+vxmejYTjdRPcBI4fT++9ezDGVMVKsuflUbx6NYXffueOYSgUHRKDSc8t/xxH8llRAGQcrJtJsKPhnYpAuD7sPuPObbxSK1B7RtJo8FKjr5ZG84kTML2e8BXV+enJ03aYuPhzIv/8cK27LYvVrlCcCZxaDVqz+AB5J9wbrbS18UpFUF9CF1nrRRgUEcmlj8xh8u33tpZYDRLbqy+XP/Fs4xVrYzTDkOvgyRy49G249Re46mN4KAU6OTep01bD0kfrVQaG0FA63XILPTZuIOoxLduoLS+vJUNRKM4IqmcxW/D0OlJ3ZDdQu33jlYqgvlhDhbXyDV815wWShpxVafLZ1gghSBgwyPXQFrXRG2DgVRA3DHpfCIGd4ZxqS17r5sEzIdo+Qml+3eYBAZUB7Yp+WErhsmXNk0OhOEOISQ7hyserHEy/m7u9zgdlR8ErFUF9ISYyD9U0Fw0Mb3kCnHZP/yugWjyTSl5MgJ+f1kJiV8PQqSqX6vF77sVu6djTYYWipUTEB3LerVUmpQUnO2aiKK9UBPUFnQuNrprmTbv3zx43B2sXCAFTX4QL/wN+nWqWrfqPdlQj8Pzz8BtVldv10AUXUH7ocGtIqlC0W7oNiaw8/3jO2lbzLXAnXqkI6psRTJ41mxGXXsmd//2Y3mPOaQOpmsb+9Wsoys1GOhxIKSmzWBpvdDqG3Qx/OVTX7HTlP+H5LuCMjSSEIPLBhyqLbZmZHJo2jbKUlOY/W6Ho4Oh0gq4Dq2bLJzqgFZF3OpTV40cQFB7J2KuvbwNpmobdGVL6+1dfqlPWbdgIxlx5LREJXZvXeewQeLpA8zX4+i7tXnkh/PQEDJ8FwXGY+/cjbu7rHLt7dmWzwxdfQu+9rZ9wW6FoL0y7cwCLnltP9lELu1enE9srFH0TzdTbko4jqRtpahjq9kRDuRQOblzHqkUftfwhg2bCiDurrtfNg9eGaD4H2QcInDiRXrt2Yh46tOXPUijOEM6Z0ROdTnBwcxbz7l5BqaWirUVymY77RmwBTU1M056YfNtsuvQdgH9oWL3lxXl5WHJz2LNqBTZrA05nDSEETH0B+tbj0v76UHg6GPG3MBKSvqmqf/h3bQnJ4YCsFNj0AZTmgcOu3T+2UbuvUJyhdE4KZtRl3Sqv33t4FXZrxwjt7pVLQx15RtDv3Mn0O3cyoEUnFTod5SXFrPzkQ7b99B3FeTnMf/guyouLmSplyxziLnsbpr8KxVnw6uC65ad+jNIBH1xUt/yb0/hgPJIKZs9HbVQoWptBk+JJ35/P4W2aT0H2MQtRXYPaWKrG6bhvxBbQkWcE1Tml0Hz8/Ol+lmbNY8nLpbxYM+ssLSxomV2z3gg+gRCWpEU7veiVFssMwIuJmr/CN/fVGw5boejITLtzAP3P0awQ139zqI2lcQ2vnBE0J3l9e+dUCAqT2UxAWDi5x4+y4sN3KMw6yfgbbsNaVorJ3IJ8CkLA0Bu14xTp22DR1SB0cNX/tOWfPtMhvAd8/2foMhzKiyCyD4QnwysDa/a5ab52AHTqruVWSDoX/GuZsioUHYzO3YLZ8dtxjuzOZceKY/QfH9fWIjWIVyqCpkYf7QhEdu3GlXOeJzQ6lnVfLiL3+FEANv+whM0/LAHgT489Q+IgN27wdu6v/SslWb8cQdrDKHnnJQwx0ZRu3kvE/VMIufTmqvpP5sCn10NKPUHrcg5owfDqo9sEOPirpmCumA9RfaE4G3Z8pt23ZEKGM+tarwth77dwy8/QpelpRRUKd9Cld9Ue3uaf0uh3Tmy79k3ySkVwJs4IhBB06aO9mPuePZGUP1ZRWljTnjnn+DGiuiWTfSSV7GNH6DNuAhn791JWbKHX6LNb9Pzs16sFtnOmwsx47DH8h5+FwRnBVOgNMGOBVsfhgD1fw2//gJO7G+78oDMjW/Y+eHN0w3X3fqv9++6kqnsT58C4B10dikLRYsyBJm76x1je/8sqLLnlZBwsIKZ7A3lE2hjR0WJjDBs2TG7cuLHJ7X5+9022/aR9id759v/wC26/vxR3UGop4rNnHiPrSKpL9W9/8wMCwpq2JCOl5MgNN1Kyvp5w107MgwZRum0bhvBwui75GkNoA5vElizYvkjzWzj1ZQ/QeQDkHoKKepzm/CO1RD3Z+5wPDIPS3IYFv/oT6DWtkdEpFC3n27nbSNuRg8Gk4/ZXx7epLEKITVLKYfWVeeWMoCNbDbmKOSCQ61/SvtJXL/qItV8sarC+taK8yc8QQpDw4Qf1lh04/3ysaUco3boVAFtWFvtHjcYQEUHCRx9iSkys2yggAkbP1o76SPsDtnwEQ67Xlod8TpMsqDQf9i2FwnT45Zm65QtnVJ2f/ReY8EQDo3QiZc20oAqFC8R0DyFtRw62CgeHtmSRNLh9xjDzyhnB7PcX4ePn727R2i3F+XnsWbmcjP0p9B0/CZ3BwN5Vv9Glb39WLphPcX4ecX36kdBvEElDhxNZLUdycznx7HMULVtGwLnjyV9YUwmZhw1FILCmpyOMRsyDBhH99+fqDQ/eYhx2zYz1+ObGl6BOMehacFhBb4JDK6DgaFVZSAIExcKRNdp1/yth0tMQHFtPRwpvx1Zh54PH1lBWbEUImH7/YOJ6to3pdEMzAq9UBPd88BkmX7O7ReuQvHHbNXX2EibcfAeDzrvAbZtbjooKsl97jZz/Npxis/NfnyH0yivd8szTknsI1r4JG98Dh5uCgw2+Dqa/pmYMinqxWx38/MFuDmw8ycCJXRh7RXKbyKGWhmpRX6whbyVpyFnsWvFzjXu/vjeP7LRUxlx9HX5BwS1+hs5kIvKhhzAldiXjiaplGP9x4yheubLy+sRTczjx1Byt7Oxx6AODQKdD6PXY8/LodPss/IYMaZkwYUkw7SXtACjJhX0/wtG1YCuHbZ+A0EOPKZA4BvwjtOPYBtj+KcQOhYgekHMQtn6s9bHlI+0Y9xBMfKpl8inOOPRGHVGJQRzYeLLdRib1yhnB/R9/1fwEL2cYUkqEEOSmH2fR049QUlCVlGbsjBsYcckVbn2eLS8PfUhI5WzDnp9Pwbffkfmsa9nXAsaPJ+711xDt4fd3dD28O/n05VNegL6XakmAFF7N1p+PsPrzAwBc/eRwOsUGtLoMakYANdIwnonmo83l1As5LCaW2+d9wI9vvMzulVo4aluF+4Nm1bYa0oeEEDpzBvrQEMpT9mHNyKDwm29q1PFJ7k75fu0/kWXFCiyrVhE4frzbZWsyXYbDXw5r/g+nTFyrs/RR7TiFTzDYSsHgC/YKCInX8kBE9ob4Udr+Q8wgMPiozekzjLheVX4FC/+2nlGXdmPI+QltKFFNPDojEEJMAV4B9MA7UsoXapVfAzzivLQAd0optzXUZ3NnBMv++zrbf14KwEOLvm1ye29BOhys+XwBaxcvZOSfrmb0Fde0C0eYspQUDl98SeV10LRpBE6eRMD48ejM7WC/x+HQnOIytsEXt7q376TxMMppSSUEdJuolEQHxJJXzgePra68vuaZkYREtcDbv4m0yYxACKEH5gKTgWPABiHEEillddONw8A5Uso8IcRU4G1ghCfkkSqmjUsInQ4htBnT2sULKcw6ydS7294Zy7dnT0KuuIL8zz4DoPD77yn8/nv8R48m/r13kTYbjrJy9AFtZA2m02l7BxE9YIBzOa00Dza8q52X5Gj+EFl7NSuknYs1qyS7C7OuQyu0ozZ+nbS9jKi+2rnRD3o7g/8pRdHuCAj1qXQyA/jy35u56cWxbSyVhieXhoYDB6SUhwCEEAuBi4FKRSClXFOt/lrAYwE5HHalCFzFYDJVnu/+/Vd0ej2Tb5vd5sH6wu++C1tuLpZffqm8V7xmDXt69a681oeFEfvPlzDGxWGMjdU2mvPzKdm6FcsvvxJw7ngCJ0xoHYHNoXD2w/WXXf5e/fctJzU/iGMbtA3slKVQfBICY6AovWbdkpyqDet6EdBzmmbeGhSjLUnpvWc1uD3iF2Ri6JQENi1No6SggtyMYsKi296U3WNLQ0KIy4EpUspbndfXASOklPV6CwkhHgZ6napfq2wWMAsgPj5+aFpaWpPl+eH1f1WufauloYYptRSxeuGHbFv2Q+W9yMRuXHDfnwmLaR/BsxwlJaQMaV7cpM5/+yuBkyY17OXcXrFkwXcPaBvVlkztnqszC4DJf9WUQ3jbmDAqwGa189Y9v1VeX/fsKILCPb+82SZ+BEKIK4DzaymC4VLKe+qpey7wBjBWSpnTUL/N3SP47tWX2Lta++ErRdA4DoedpW+8zJ6VVXmMh110GVFdu7H1p+8pLSzAYbeTn5mBwceHIVMuYtCUCwkMC2+gV/ditxRTsn4djuISfPv24fj9D1C+b1/digaDtoZfa3kw4ZMF+A2uJ89CR0RKyEvVvK1L87XZxNd3u9ZW6OD85yE4DnR60Bm1/Q6HVZtJdBmpzSb8wrTcE0Knlp5ayP4Nmfz07i4ATGYDt/2nZbG+XKGtrIaOAV2qXccB6bUrCSEGAO8AUxtTAi3BofYImoROp2fa7IeI7t6DX99/C4CN33xRb11beTnrv/4ch8PBOdfeXG8dT6AP8K+xzJO05Gukw4EtKwvrsWMYY2KQNhvGqCiEyUT2vLfIevnlyvppM2aSvGY1hrD6s711KISAMGeuav9wCO8Og6/VrsstsOZV+O3F+ttKByx9pP6yhvAJhnKnM+LYB5zLVxmQMAaCoqFTMhhMDffhpSSfFUXGgXx2/HacilIbR3blEN+37cKve3JGYAD2AROB48AGYKaUcle1OvHAr8D1tfYLTktzZwRL/v139q/THqFmBE1j8w9LWD7/7Tr3TWYz0iGxlpcBMHjKRUy46fbWFq/JZM2dS/ZrWhwmXUAASd99i/XYMYrXrsXvrLPwHz68jSX0IFJq+w/WEqgogUXXaIrgFKFdtRAaaavc87wx94NvsOaIF9lH26NQ2emAmuEnfAOMXPPMSHz9jR57XpvMCKSUNiHEbOBHNPPR96SUu4QQdzjL5wFPAZ2AN5wmirbTCdpS1GZx8+l37mQMRhPRPXoREZ9Yp/yUokj5YyUlBfmUFVsYfcU1BEdG4RsQ2O6c98LvuovCJd9QkZaGw2LhwDnj61bS6ej81JPY8/Ox5xfQadZtGMLCkHY7FWlpVKSl4duzJ/bCQmRFBb79+nWMYIZCaP4Pp5iT13ibUx+L5YXaslHmLjCaITBaC+q3+QPoPklbPtpb6yNr9ct1+ztF5/5arKb0LRCaCLkHtT5Du2oJkMoKNGuoM3SD22DSc9M/xrDouQ3kphfzx5cHOffaXm0ii9d4Fn/54jMc2rwBUDMCd3O6GQNATM8+zPjrP1pZosZxlJdzaOo0rOnp2h6CzUXX/wbqCh8foh59hNAZM+ot9yp+fQ7SN8OJnWA50bK+DL5gK9NmFJ26wZ5vYOhN0HOqFjLEaNaWpTqCIq6Hr/6zmeMp+ZiDTNz8D8+Zk6qgc8Di5+eQunUToBSBu8k8fJAvnp9DSFQ06fv31vDiBhgwcQpZaYcZc9V1JAwY1EZS1o+jpARhNCKMRmy5uZSnpJD/xZcUfvMNfsOHU7Z3L47CwprlXVFxAAAWPElEQVSNjEawWkGvB7u9wf714eH49u5N8PTpBF90oQdH0kEozYODy7WXucFHyzQXmgAn92rlp5akhK7mkpWr9LpQ89cYe7/WfwfgxKECFv9Dezdd9vAQoj2UwEYpAuDz554kbfsWQCkCT2OtKCc/I50P/1LTQKzfuZM5/4772kiq5lNx5Ahlu3ahCwjEPKA/+uBgpNUKQmAvKiJvwQLK96ZQtGxZg/2E3zObiLtdtOTxZk4ZdthKIXWV5kxnyYTogZrZ7P5lYHchf0ZkXy1kx6HftGWouKFaKA9TgNZXO7F8slsdzLtnReX1rFfPwWhyv8+OUgTAp399nKO7tgNKEbQGUkqWzv03+SczMfn6krptMwA9Ro3jovubYaHSAXCUlFB+8CDY7ThKS8maO5fylH04iooq6wScey5xr73aPoLmnSlICev/qzncrfpP09pevUDzq2hjpVBSWFHpcQxw9zz3Oz2qoHOoEBOtjRCCqbMfAiB166ZKRbB/3eqGmnVodH5+mPv3r7z2HzUKgNLt20m98ioALMuXs7efVkeYTEQ8+ADG2FgCJ2rxg9pDXKcOhxAwYpZ2Pulpbelpx+eaOW36FshLg8wdpyoD1T5+F86sOo8dps084kdqXtsX/Evbg2gF/IJMxPYM4XiKFv33wKaTdB8a2SrPBi+aEXzy5J9J37cHUDOC1kZKSdqOrSx+7kkAgiKimPnsP/EP8R4zQltuLgcmTUaWlDRa1yc5mfL9+wEInz2b8Fm3YS8qAp2uRghvRTNx2LXcEl/d0XjdoDjNsS7uLC03dmRfLbZTcJzmr2Er1/JVuOl3MveOqii2Vz81nE4x7gtXrZaGgAVPPETGgRRAKYK2wG6zMe+O6ykr0jZeEwYMJrp7D7KOpHJw4zqGXXQZAydNJSgyEofdgcHoOXvqtsRRUkLZ3hRy33+f4rVraywbtYTQa67B1C0Jn8RE/IYNQ5iUI5dL2Cq0ZEQ6vTZz0BvhyFo4+EvjbasT0RvOugVihmjhxX0CobxI27D2CXRZUaRuz+a7N7Ql7GHTEhkx3X0zEqUIgP89dj+Zh7SY9koRtA35mSf4+PEHKLOc/uVnMJqQSK574RU6xcVjt9mwWyswmVsvXG9r4ygvBykp37sXa1YWxugYyvfvJ+Oxx5rdp3nIECIfehBDeDimhPYT977DYC2FA79A9j4ttLjDpvlI+AZr/g3NIX4UZO+HiJ4w4g7NA9u/rjfx4W1ZfP+mtpQ15fZ+JA2KcMssUCkC4MNH7iUr9RCgFEFbsnP5Mn776F3Kii30n3g+SMmOX386bX2dXo/Dbsc3MAgBBEdG0WvMeMJiYonu0Qujjy/5J9LJzzxBdHJPt6TWbE/YCwvR+fsj9HocxcVYM0+C3YYtJxfhY6Lop2Xkf/opvn37UrJ+fb19JC7+HN8ePRBn6CyrTSlMhx+fgP0/QYWl6n5TAgGeottEmPI86cd1fDnvaOXtqXf0J2lQRItFVYoA+ODhu8k+qkUtVYqgfSGl5PDWjRiMJnb8+lNlcMDm0OfsCUy58/6O4eXrAWzZ2aRdfwMVhw41XlkIkBK/USOx5xdgTUtDFxxMzAsvYEpMwBAejmjj0OMdkurZ5exWOLEdcg9rx6b3ofB4g80dUseaohvYVjIdgDFdVzNoUBkMuV7ztPZvXmBHpQiA9x+8k9zjmpZViqD9Umop4uiu7RSczCQyMQlf/wD2r/+D/MwMso+kknPsiEv9xPcbyNnX3ERUUncPS9x+Ofnv/5Dzzjt1oq42hR7r16EPCnKjVIpKLFmQsVWLGvt93bwV24ovYFXRrRhFKTPC7yFQn6PFbLqtnrSoLqAUAfDe/bPIy9CCnypF0PEpzDoJAsxBwRhNPhzbvZPFL8zBVl7T0cg/NIweI8eQe/wYxXm59J84hf4TJlNRWgpoSXiEEGfsHoR0OEAIZHk5jqIibLl56PzM2E6cwG6xYM/NxRARgT44mJx33qFo2c91+gi99lrKU1LQBQTgKC7GPHgwsqwMvxEj8OneDX1YJ+z5+WCzYkxIUFZN7sCShX3zxyxe2p2sXO1v8/Zu92GI7QMzFzWrS6UIgHfuvZWCTC3miVIEZyZSSrLSDvPbR+9yZGeDqa9PS3T3nvQaew5hsV2I7t4DH7/Gs0dJh4NSSxHmgMAzZknqyKxZFP++slltgy68EF2AP0HTpuGTnNwxEwC1E7KPFbHoWS1G2tlX96D/+OYnhlIOZSiHMm9ACEFkYhJXPPkc+9et4cDGtRhMJkI6x3B87y4OblzXaB8ZB1IqzYwBwmK7EBQRicnHF6HTYcnLJfPwAYTQYS0rJTA8AktuTuXf16jLZ+IfEoolL4fO3XoQ338gRlPHiHlTnZi//53sN95AGI1YT2Tik5yMZeXv6Mx+2LKyED4myndrfjnCbEY6Z1gAhd9qH1r5C2t+ufqfPY6wG26o9IUQPj6YunQB50a4o6gI4eenFEc1wuMC0Rt12K0OjqfktUgRNITXzAjeuvMGLLla3hs1I/BOpNRyJ2jLQdqXu8Nuozg/jyM7trHr9184efhg5bKRu7j8iWfx8fMjJDoGX/8ApJSUFVvw8fNDp9MjpeyQyylSSmRFBTofH6TVSsGSJZTt2k3p1q1Im63+bHEuEjT9Isz9+mHLzQW7HeHri0+3bvh07165lCWl1MrO8HAdp8xJgyLMzHhyOIZmxiFSS0PAvNuvozhfi72uFIGiMfasWsHBjesw+popL7EQ3b0nAaFhmIOCyU0/hsFoIiyuC4FhnfALCWXdF5+y9adv6RSXQGh0DLtW1F1rrw+T2Y+K0hJMZjNJQ4YTldQdS24OwZFRdD9rFIGdWi/1p7uRDofmQLdzJyeefoaK1FT0YWHYc3PRd+qEPacqIaHO3x9HcbHrnQuhhZ12Rn81Dx5M2Z49GCIjEUYjjuJibCdO4NOnN9a0I9rexsCBWNPTiXricYKmTHH3cD2G3epg0d83kJdRzODJ8Yz+U/MMIJQiAN64dSalTq9WpQgUnkZKybZlP7Dm0/9h9PXVNrebgW9gUKU3tm9gEP3GT6Jzt2R8/QMJ6dwZv6AQjL6+7hS91ZAOB9Jq1cKAO/dWLCtXkvWflzF164bO34/Sbdu1fzduQvj4IMtdiDrqIvrQUAInTSL8zjswxsS4rV9PcOJwAV/8YxPBkX5c9X9nYTA2fVagFAEw9+arKSvWHD6UIlC0NnabjTJLEXarFd/AQHR6A5bcHEy+vuxc8TN7Vi4nLC4ek6+ZgxvXVn60uIrJbKb/xCmUFRVhycshOCKK6OSeAPQcNa7DKovaSCmxncxCZ/ZFZzZTtncv9oJChNGIrKhAHxyEzmzGXlgIUiKMRipSU5FWK7bsnBo5qxvCd8AAZGmp5k/RORqAgLFjMA8dij7AffF/msLhbVnE9Q5rdohqpQiA1268onLtVykCRXvHWl5GfuYJfAMCKMzKYv+6VaTt2IYlN6fBEB2no+vgYWQeOkB8v4H4h4TS5+wJRCa2TmTN9oSjrIzitWsp3byF/MWLayxPuYrfsGE4yssJOv88DBEROEpLMSUloTOZ0AUHow8IQPj6og8M9MAImo9SBMAr1/+p0sZcKQLFmYCUkvLiYrb++C3lpSWYA4Pw9Q9g03dfEdI5ujI1a2P0HT+JyISulJeUkJ+ZgcnsR5e+/emclIzBxwdf/wB0Z7CHsaO8HEdJCbLCSvm+fQijgdKt25AV5ViPp1O+fz9lu3c3vWODgeDp06lITcV/5Eik3Y49Px+fHsnICiv+Y0bj0717q5kcK0UAvHzNJdiduWaVIlB4AxWlJexf/wcA+SfSKcrJpig3hyM7tjarP9+AQEKjY3DYHRTn5dBt2Ei6Dx+FtbwMa1kZ5SXFRMQnYrfasFkrKC0qxC84GGtZOX7BwYTHJ+LrH9AhLaQArBkZlO3ZQ8E331CRmoZPUhKl27ah8/enPCUFQ0w0tvSMJvfrP2YMgeefp52PHoMpLtbdogNKEQDw7xnTK229lSJQeDO2igoqykrJTT/GlqXfYjAa8fHzx+FwYMnNdsnfwh34BYcQEhWNdDgIi+1CbsYxgjpFIHQ6jD4+hHdJwFpRQXT3HvgFBWP09SU4snOryNYSpNVK0c8/Yy8soiI1FRx2hJ8fZdu2YYyPp+JwKiXrTv8z9u3bF2mzIe02Kg4cxKdPb2RFBVGPPYZP167N3thWigD411VVicOVIlAoXMNht5Obfozso2kYTD74+PlxdNd2Dm5aj91q1V7QZjOZB/fj4x+Ab0AgBqOR7KNpBHaKwBwUROrWTW6VSW8wICVc/dcX6dytR2WQt44206hISyPz78+jCw5ClpY1mvMatPDiiQs+btbzlGdxNU45EikUisbR6fWEd0kgvEtVToMuffoz+oprmtyXzWrFYbdRmHWSrCOplBYWUlKQR2CnCMosRZSXFJNz/CgVpSWERsfisNnY9dsvBEdGUXAys7KfU0u8C554qJasBkDisNvp3C2Zwuws/ENCCYqIIvPwAbr07ofNWkH+iQx8/QOI7d0Xa3k5+ScynEtedhx2Gw67nYKTmXQdNBS7zUZ8v4F07pbcvB9gA5gSEujy1rzKa2tmJtb0dITBiDDokXYH9rw8LMuXk7dgAcaYGAyd6uYvcAdeNyPQ6Q08sOArd4ulUCg8TKnTWmrn8mX8/r/32lSWiISuWMvLKC0sJCyuCw5blRLJOX6UHiPH4rDZcNhtZB89QkyPXtitVkoKCwiLjaMoJ5vATuGYfM2UFBYQ33cAySPGYDKbPTazUUtDVCkCg8mH+z5a7G6xFApFK1NRVopAYDCZsFaUIxAUZGVSVlSERJJz9Ag+/v6UFhagMxgx+viQe/wopUWFmMxmTGZ/8jMzCIuORWcwoNPr0ekNZKUdQjokBpOJbcu+b9Ux6fQGHHZtxpMwYDBhsXGUWSwk9B9EdHJPDCYTQeHNS2qvloaqcaZEh1QovB2Tr7nOee0lrJYy8eY7sFkrKLNYKCu2YLda0em1+FAl+XkYTCZ0egM6g6ZECk6ewFpWht5gQGcwUJKfj8nPD6SkICsTv6AQSosKkdKhhU/fu4vMwwcpyc9DCFG57AWQtn0Ladu3ALBn5XIAonv0Yubf/tnicdXG6xSBTikChULhIpoFky9GH1+X4j5Fde3WpP6HXXRZjWtreRkFmScqw6in7dhKVlqq0xQ3hICQsCb17ypepwhU6j2FQtFeMfr4Eh6fSHh8IgBDpl3cKs/1us9jNSNQKBSKmnjdW1EpAoVCoaiJ170VhU4tDSkUCkV1PKoIhBBThBApQogDQohH6ykXQohXneXbhRBDPCkPgE7vdbpPoVAoGsRjb0UhhB6YC0wF+gAzhBB9alWbCiQ7j1nAm56Sp1IutTSkUCgUNfDkW3E4cEBKeUhKWQEsBGpvgV8MfCg11gIhQohoD8qEsZrtsUKhUCg8qwhigaPVro857zW1DkKIWUKIjUKIjVlZWc0S5rLHniEioSsX3PNws9orFArFmYon/QjqC5hRO56FK3WQUr4NvA1aiInmCNN10FC6DhranKYKhUJxRuPJGcExoEu16zggvRl1FAqFQuFBPKkINgDJQoiuQggTcDWwpFadJcD1TuuhkUCBlLLpKX4UCoVC0Ww8tjQkpbQJIWYDPwJ64D0p5S4hxB3O8nnA98A04ABQAtzkKXkUCoVCUT8ejTUkpfwe7WVf/d68aucSuNuTMigUCoWiYZRRvUKhUHg5ShEoFAqFl6MUgUKhUHg5ShEoFAqFl9PhchYLIbKAtGY2Dwey3ShOR0CN2TtQY/YOWjLmBCllRH0FHU4RtAQhxMbTJW8+U1Fj9g7UmL0DT41ZLQ0pFAqFl6MUgUKhUHg53qYI3m5rAdoANWbvQI3ZO/DImL1qj0ChUCgUdfG2GYFCoVAoaqEUgUKhUHg5Z6QiEEJMEUKkCCEOCCEeradcCCFedZZvF0IMaQs53YkLY77GOdbtQog1QoiBbSGnO2lszNXqnSWEsAshLm9N+TyBK2MWQowXQmwVQuwSQvzW2jK6Gxf+toOFEN8IIbY5x9yhoxgLId4TQpwUQuw8Tbn7319SyjPqQAt5fRBIAkzANqBPrTrTgB/QMqSNBNa1tdytMObRQKjzfKo3jLlavV/RouBe3tZyt8LvOQTYDcQ7ryPbWu5WGPPjwIvO8wggFzC1tewtGPPZwBBg52nK3f7+OhNnBMOBA1LKQ1LKCmAhcHGtOhcDH0qNtUCIECK6tQV1I42OWUq5RkqZ57xci5YNriPjyu8Z4B5gMXCyNYXzEK6MeSbwhZTyCICUsqOP25UxSyBQCCGAADRFYGtdMd2HlPJ3tDGcDre/v85ERRALHK12fcx5r6l1OhJNHc8taF8UHZlGxyyEiAUuBeZxZuDK77kHECqEWCGE2CSEuL7VpPMMroz5daA3WprbHcB9UkpH64jXJrj9/eXRxDRthKjnXm0bWVfqdCRcHo8Q4lw0RTDWoxJ5HlfG/DLwiJTSrn0sdnhcGbMBGApMBMzAH0KItVLKfZ4WzkO4Mubzga3ABKAbsEwIsVJKWehp4doIt7+/zkRFcAzoUu06Du1Loal1OhIujUcIMQB4B5gqpcxpJdk8hStjHgYsdCqBcGCaEMImpfyqdUR0O67+bWdLKYuBYiHE78BAoKMqAlfGfBPwgtQW0A8IIQ4DvYD1rSNiq+P299eZuDS0AUgWQnQVQpiAq4ElteosAa537r6PBAqklBmtLagbaXTMQoh44Avgug78dVidRscspewqpUyUUiYCnwN3dWAlAK79bX8NjBNCGIQQfsAIYE8ry+lOXBnzEbQZEEKIKKAncKhVpWxd3P7+OuNmBFJKmxBiNvAjmsXBe1LKXUKIO5zl89AsSKYBB4AStC+KDouLY34K6AS84fxCtskOHLnRxTGfUbgyZinlHiHEUmA74ADekVLWa4bYEXDx9/w3YL4QYgfasskjUsoOG55aCPEJMB4IF0IcA+YARvDc+0uFmFAoFAov50xcGlIoFApFE1CKQKFQKLwcpQgUCoXCy1GKQKFQKLwcpQgUCoXCy1GKQKFoRYQQiaeiSjqjhH7b1jIpFEoRKBQu4HTeUf9fFGck6g9boTgNzq/3PUKIN4DNwHVCiD+EEJuFEJ8JIQKc9c5y5njYJoRYL4QIdLZd6ay7WQgxum1Ho1CcHqUIFIqG6Ql8CExGC9Y3SUo5BNgIPOgMe7AILeLlQGASUIoW9nqys+5VwKttIbxC4QpnXIgJhcLNpEkp1wohLgT6AKudITpMwB9oiiJDSrkB4FTESyGEP/C6EGIQYEcLD61QtEuUIlAoGqbY+a8AlkkpZ1QvdEZ0rS9OywNAJlrkTx1Q5kkhFYqWoJaGFArXWAuMEUJ0BxBC+AkhegB7gRghxFnO+4FCCAMQjDZTcADXoQVMUyjaJUoRKBQuIKXMAm4EPhFCbEdTDL2c6ROvAl4TQmwDlgG+wBvADUKItWjLQsX1dqxQtANU9FGFQqHwctSMQKFQKLwcpQgUCoXCy1GKQKFQKLwcpQgUCoXCy1GKQKFQKLwcpQgUCoXCy1GKQKFQKLyc/wfdUNujqOZhSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# precision recall curve\n",
    "precision = dict()\n",
    "recall = dict()\n",
    "for i in range(6):\n",
    "    precision[i], recall[i], _ = sklearn.metrics.precision_recall_curve(y_test[:, i],\n",
    "                                                        y_score[:, i])\n",
    "    plt.plot(recall[i], precision[i], lw=2, label='class {}'.format(i))\n",
    "\n",
    "plt.xlabel(\"recall\")\n",
    "plt.ylabel(\"precision\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.title(\"precision vs. recall curve\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
